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1.        INTRODUCTION 

In  recent  years  a  study  of  the  gasdynamics  of  space-based  HF  DF  laser  exhaust  plume  was 
conducted  at  the  NPS.  The  aim  of  this  program  was  to  identify  phenomena  that  give  rise  to  a 
potentially  contaminating  backflow  of  corrosive  species  (HF,  DF,  F),  and  to  estimate  the  magnitude 
of  the  flux  arriving  at  the  spacecraft  [1,  2,  3,  4]  .  The  laser  was  envisioned  as  a  cylindrical  spacecraft 
having  a  centrally  located  nozzle  of  ring-symmetry.  This  is  an  idealization  of  a  more  general  zero- 
thrust  exhaust  configuration  for  an  open  loop  HF/DF  laser  (Fig.  1-1). 

In  former  studies  [2,  3]  we  focused  on  the  contribution  of  thermal  self-scattering  to  a 
contaminating  backscattered  flux  of  corrosive  molecules  (HF,  DF)  from  the  exhaust  plume.  It  was 
shown  that  this  flux  emanates  primarily  from  the  lip-centered  rarefaction  fans  that  flank  the  exhaust 
plume,  and  it  was  estimated  that  at  a  sufficiently  high  exit  Mach  number  (e.g.  M.  =4),  the  thermally 
backscattered  flux  of  HF  +  DF  is  utterly  negligible. 

The  purpose  of  this  report  is  to  present  a  study  of  the  effect  of  diffuser  start-up  flow,  including  the 
ongoing  recombination  of  atomic  hydrogen,  on  the  thermally  backscattered  flux.  Our  main 
conclusion  is  that  by  designing  for  a  sufficiently  high  steady  exit  Mach  number  (e.g.  M,  =  4),  the  level 
of  thermally  backscattered  flux  of  corrosive  molecules  (HF+DF)  would  be  negligible  even  when  the 
combined  effects  of  hydrogen  recombination  and  transient  flow  are  considered. 

The  diffuser  is  simplified  as  a  cylindrically  expanding  channel  of  uniform  width,  and  the  flow  is 
idealized  as  inviscid  expansion  of  a  perfect  gas  with  an  ongoing  reaction  of  hydrogen  recombination. 
The  flow  is  computed  as  one-dimensional  time-dependent. 

A  fully  2-D  time-dependent  computation  of  an  emerging  exhaust  plume  is  outside  the  scope  of  the 
present  laser  exhaust  study.  Rather,  we  resort  to  semi-quantitative  arguments  that  lead  to  an  upper- 
bound  estimate  of  the  effect  of  transient  flow  on  thermal  backscattering,  by  showing  that  the 
transient  turning  angle  in  an  expansive  corner  flow  can  be  estimated  as  higher  than  the  corresponding 
steady  flow  turning  angle.  It  is  subsequently  suggested  that  an  overestimate  to  the  the  backscattered 
flux  from  a  transient  plume  is  obtained  by  considering  it  as  a  steady  plume  whose  exit  velocity  vector 
is  pre-rotated  so  that  its  limiting  (vacuum)  streamline  coincides  with  that  defined  by  the  transient 
turning  angle. 


The  plan  of  this  report  is  the  following.  In  Ch.  2  we  present  our  radial  (1-D)  diffuser  model  based 
on  the  linear  configuration  of  the  TRW  test  laser  [5]  ,  and  we  propose  a  "typical  case"  of  HF,  DF 
laser  How  based  on  one  of  those  tests.  Ch.  3  is  devoted  to  the  gasdynamic  governing  equations, 
including  the  hydrogen  recombination  rate.  The  resulting  dilTuser  model  was  implemented  in  a  1-D 
Euler  code  (named  HFL)  which  utilizes  the  GRP  scheme  for  integrating  the  conservation  laws  of 
compressible  flow  [6]  .  The  HFL  code  is  given  in  Appendix  A.  The  results  of  a  typical  dilTuser  start- 
up flow  are  presented  in  Ch.  4,  along  with  a  discussion  and  analysis  of  the  effects  of  transient  flow 
and  uncertainty  in  the  hydrogen  recombination  rate  on  thermally  backscattered  flux  of  HF+DF.  We 
also  consider  kinetic  backscattering  of  HF  or  DF  molecules  caused  by  third-body  recoil  in  the 
hydrogen  recombination  reaction.  It  is  shown  that  the  combined  contribution  of  these  effects  to 
backscattered  flux  remains  negligible  in  the  typical  case.  Chapter  4  ends  with  some  concluding 
remarks,  and  references  are  brought  in  Ch.  5. 
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2.        THE   DIFFUSER   FLOW   MODEL 

Our  diffuser  model  is  based  on  a  schematic  radial  laser  configuration  as  shown  in  Fig.  2-1.  The 
major  components  in  this  design  concept  are  those  present  in  an  experimental  HF/DF  open  loop  laser 
tested  at  TRW  [5]  .  In  this  chapter  we  describe  the  radial  configuration  including  the  diffuser 
(Section  2.1),  then  we  present  a  typical  HF/DF  laser  flow  based  on  one  of  those  TRW  tests 
(Section  2.2). 


2.1      Radial   Laser  Configuration 

We  assume  that  an  open  loop  HF/DF  laser  of  the  type  tested  at  TRW  [5]  ,  can  be  rearranged  in  a 
radial  (rather  than  linear)  sequence,  where  the  flow  begins  at  the  hub  and  proceeds  in  an  outward 
radial  direction.    Referring  to  Fig.  2-1,  the  major  components  of  this  configuration  are  : 

(a)  Combustion  Chamber  :  Deuterium  and  fluorine  burn  with  excess  fluorine,  resulting  in  hot 
gaseous  mixture  where  fluorine  is  virtually  completely  dissociated. 

(b)  Nozzle  Cascade  :  Rapid  expansion  to  supersonic  flow  leaves  atomic  fluorine  concentration 
effectively  frozen. 

(c)  H7  +  He  Injection:  Mixture  is  injected  between  the  supersonic  streams  of  DF+F  emerging 
from  the  cascade. 

(d)  Mixing  and  Lasing  :  The  lasing  is  from  vibrationally  excited  HF  molecules  produced  by  direct 
reaction  between  H^  and  F.   As  a  by-product,  one  H  atom  is  produced  for  every  HF  molecule. 

(e)  Diffuser  Entrance  :  This  point  marks  the  end  of  the  lasing  process.  From  this  point  on  the 
flow  is  just  an  exhaust  to  be  discarded  safely. 

(0  Radial  Diffuser  :  The  purpose  of  the  diffuser  is  to  raise  the  flow  Mach  number  at  the  exit  so 
that  no  appreciable  backflow  from  the  exhaust  plume  will  take  place  [2,  3,  4]  .  It  should  be 
noted  that  a  desirable  diffuser  area  ratio  can  be  achieved  at  lower  radius  ratios  by  letting  the 
diffuser  expand  in  the  axial  direction. 


t  -> 


Typical    Laser   Flow 


We  chose  as  a  typical  laser  flow  one  of  the  tests  conducted  in  the  TRW  experiments  [5]  it  is  test 
III  in  Table  5  (p.  91)  of  that  reference.  The  mole  fraction  and  corresponding  average  molecular 
weight  and  specific-heat  ratio  (W  and  y)  at  the  diffuser  entrance  for  this  test  were  : 


Mole  fractions 


Average  molecular  weight 


fH     =    '°91        fHF   =    '°91 


f  =.104 


fDF  =    -135       fHe 


VV  =  7.27 


.579 


(2.2-1) 


Specific-heat  ratio 


y  =   1.54 


The  additional  data  required  for  the  computational  model  is  the  flow  variables  at  the  diffuser 
entrance.  Some  of  these  variables  have  been  measured  directly  or  indirectly  and  the  remaining 
variables  can  be  evaluated  from  standard  isentropic  flow  relations  for  compressible  flow  of  an  ideal 
gas.   The  data  given  in  the  TRW  tests  report  [5]   are  : 

Diffuser  entrance  cross-section  is  :      1"  x  7" 


-i 


Mass  flow  rate  is  :       14.99      (gr  sec    ) 


(2.2-2) 


Stagnation  temperature  is  :       1400    (K) 

Flow  velocity  is  :       2300    (msec"1) 

Flow  density  is  now  obtained  as  the  ratio  of  the  specific  mass  flow  rate  and  the  velocity.  Then 
Mach  number  is  extracted  from  the  standard  relations  between  Mach  number,  velocity  and  stagnation 
temperature.   The  values  of  these  variables  are  : 


M2  =  2.26 


p2  =   1.44  x  10"3    {kg  rn3) 


(2.2-3) 


P7  =  9.72  x  10"4  (MPa) 


The  flow  at  the  diffuser  entrance  is  thus  fully  specified  for  the  selected  typical  case.  In  the  next 
chapter  we  take  up  the  matter  of  computing  the  transient  expansion  following  an  abrupt  start-up  of 
the  diffuser  inflow. 


DIFFUSER  EXIT    X  =  2.5  (m) 


DIFFUSER 


DIFFUSER  ENTRANCE 
X  =  0.625  (m) 

LASING  SECTION 

NOZZLE 
CASCADE 


PREHEAT 

D2  +F2->DFtF 


Figure   2-1.  Radial   Configuration   of  HF/DF   Open-Loop   Laser 


3.        THE   GASDYNAMIC    FLOW   MODEL 

The  transient  expansion  of  the  lasing  products  following  an  abrupt  start-up  of  inflow  at  the 
diffuser  entrance,  is  idealized  as  an  inviscid  compressible  flow  of  a  mixture  of  perfect  gases  with  one 
ongoing  chemical  reaction  -  that  of  hydrogen  recombination.  The  computational  model  thus  calls 
for  a  formulation  of  the  governing  equations  and  for  a  numerical  scheme  capable  of  integrating  them 
in  time  and  space. 

In  this  chapter  we  describe  the  governing  equations  for  the  diffuser  flow  (Section  3.1),  following  by 
a  simple  approximation  chosen  for  the  hydrogen  recombination  rate  as  a  three  body  reaction  (Section 
3.2).  The  numerical  scheme  employed  for  the  solution  is  of  the  Generalized  Riemann  Problem  (GRP) 
type  [6]  .  This  scheme  has  been  implemented  in  a  code  which  was  adapted  to  the  diffuser  flow  case. 
The  code  (named  HFL)  and  some  brief  description  of  features  introduced  to  treat  hydrogen 
recombination  are  given  in  Appendix  A.  In  a  recent  report  [7]  a  very  similar  version  of  this  code 
(without  chemical  reaction)  was  described  in  considerable  detail. 


3.1      The   Governing   Equations 

The  expansion  of  lasing  products  through  the  diffuser  is  governed  by  the  Euler  equations  for  a 
gaseous  mixture  of  ideal  gases  with  ongoing  hydrogen  recombination  reaction.  We  write  these 
equations  in  the  quasi-one-dimensional  format  of  flow  in  a  stream  tube  of  varying  cross-section  area 
A(x),  but  in  actual  computations  we  set  A(x)  =  x,  thereby  reducing  the  equations  to  the  cylindrical 
case.   (The  code  HFL,  however,   can  accept  any  smooth  A(x)). 

We  chose  to  simplify  the  designation  of  atomic  hydrogen  mole  fraction  by  normalizing  it  as 
fHZ(x,t),  where  fH  is  the  mole  fraction  at  the  diffuser  entrance  (it  is  constant)  and  Z(x,t)  denotes  the 
degree  of  dissociation  (Z=  1  is  maximum  concentration  of  H  at  the  diffuser  entrance,  Z  =  0  is 
complete  recombination).  As  an  approximation,  we  assume  that  W  and  y  are  constant  throughout. 
In  the  typical  case  (2.2-1),  y  and  W  can  vary  by  as  much  as  1%  and  4%  respectively  at  full 
recombination.  Neglecting  that  variation  is  consistent  with  the  preliminary  nature  of  the  present 
diffuser  flow  analysis. 

The  governing  equations  are  the  three  standard  conservation  laws  (mass,  momentum  and  energy), 
and   an    additional    species    conservation   law   for   atomic   hydrogen.     The   conservation    laws    are 


augmented  by  an  equation  of  state  (ideal  gas)  and  by  a  rate  law  for  hydrogen  recombination.    The 
full  system  of  equations  is  : 

Mass  [Ap]t  +  [ApU]x  =0 

Species  H  [ApZ]t  +  [ApUZ]x  =  ApZ 

Momentum  [ApUjt  +  [pU\  +  APx  =0  (3.1-1) 

Energy  [AE]t  +  [A(E  +  P)U]x  =  0 

Equation  of  State  P(p,E,Z,)  =  (y-  I)  [E  -  pU2/2  -  (pf„Z/W)QH] 

•         « 

Recombination  Rate  Z  =  Z(p,P.Z) 


where  A  is  a  function  solely  of  x,  and  p,  U,  P,  E,  Z  are  all  functions  of  (x,t).    The  normalized  rate 
function  is  Z.    Its  derivation  and  explicit  expression  are  given  in  Section  3.2. 

We  note  that  the  energy  equation  contains  the  heat  released  by  hydrogen  recombination  (Q„  per 
mole  of  H)  in  an  implicit  way,  by  defining  the  total  energy  as  including  the  latent  heat  due  to 
recombination,  in  addition  to  internal  and  kinetic  energy  terms. 


3.2      Hydrogen   Recombination 

The  hydrogen  recombination  is  commonly  assumed  to  be  a  three-body  reaction  [8]  ,  with  the 
following  rate  law  : 

H  +  H  +  M  ->  H2  +  M 

[H2]  =    -k0[H]2[M]  (3.2-1) 

kn  =  7  x  io15  [(mole/cm3)"2  sec-1] 


where  by  M  we  denote  a  third  molecule,  which  in  our  case  may  be  any  molecule  in  the  gaseous 
mixture  (including  H). 

The  value  we  chose  for  k0  is  an  upper  limit  of  a  range  of  values  recommended  by  Kondratiev  [9] 
for  the  reaction      H    +    H    +    He    (from  2*  10l:>  to  7  *  10    ),  since   in  the  typical  case  helium 
constitutes  about  58%  of  the  molecules  in  the  lasing  products. 

The  data  compiled  by  Kondratiev  [9]  indicates  a  comparable  range  of  variation  with  temperature 
as  with  the  third  species  M.  We  later  introduce  uncertainty  factors  of  10  and  100,  which  in  all 
likelihood  more  than  reflect  the  uncertainty  in  the  reaction  rate  for  hydrogen  recombination. 

It  is  convenient  to  redefine  the  reaction  rate  in  terms  of  flow  density  p  and  normalized  H 
concentration  Z.   The  modified  rate  constant  kj  is  given  by  : 

Z  =   -  kj  p2  Z2 

(3.2-2) 
k,  =  2  k0  fH  /  W2 

where  the  factor  2  is  due  to  the  fact  that  the  rate  of  depletion  of  H  is  twice  the  rate  of  production  of 
Hv  In  deriving  (3.2-2)  we  made  use  of  the  relations  fH  =  [H]-,/[M]2  and  [M]=p/W,  where  index  2 
denotes  the  diffuser  entrance. 

The  heat  of  formation  released  per  mole  of  recombined  H2  is  435.783   (kJ/mole)  according  to  [10] 
(page  F-179  in  that  reference).     In  the  code  we  use  the  energy  parameter  per  unit  mass  of  the 
gaseous  mixture  QDET.  defined  as  : 

QDET  =   QH  fH  /  W 

(3.2-3) 

QH  =  435.783/  2       (KJ/mole) 


4.        RESULTS   AND   DISCUSSION 

Several  numerical  computations  of  the  transient  difTuser  flow  were  performed  using  the  code  HFL 
(Appendix  A).  The  flow  was  started  by  an  abrupt  inflow  at  the  difTuser  inlet;  the  radial  difTuser 
model  and  the  typical  HF/DF  laser  case  (Sections  2.1  and  2.2)  were  assumed.  Due  to  uncertainty  in 
hydrogen  recombination  rate,  three  rate  levels  were  assumed  :  the  nominal  level  (3.2-1),  a  tenfold 
increased  rate  and  a  hundredfold  increased  rate. 

We  are  primarily  concerned  with  the  thermally  backscattered  flux  arriving  at  the  surface  from  the 
exhaust  plume.  Our  goal  is  to  show  that  by  choosing  a  difTuser  with  a  sufficiently  high  steady  exit 
Mach  number,  the  backscattered  flux  can  be  negligibly  low  even  when  the  combined  effects  of 
hydrogen  recombination  and  transient  exit  flow  are  taken  into  account.  Put  in  other  words,  we 
present  a  simplified  analysis  of  the  transient  exhaust  flow,  which  enables  the  establishment  of  a 
reasonable  margin  on  a  frozen/steady  design  for  low  backscattered  flux.  This  analysis  is  presented  in 
Section  4.1  below. 

The  continuation  of  the  hydrogen  recombination  reaction  into  the  plume  may  give  rise  to  a 
contaminating  backflow  of  a  different  kind.  The  velocity  imparted  to  the  "third  body"  molecule  in  the 
hydrogen  recombination  process  (Eq.  3.2-1)  may  be  large  enough  to  overcome  the  radial  flow  velocity 
component,  resulting  in  molecular  backflow.  This  effect  is  discussed  in  Section  4.2,  and  it  is  shown  to 
be  potentially  capable  to  give  rise  to  an  HF/DF  deposition  rate  of  the  order  of  1  molecular 
monolayer  per  hour. 


4. 1      Thermal   Backscattering 

Since  in  steady  flow  thermally  backscattered  flux  is  linked  primarily  to  exit  Mach  number,  we  focus 
our  attention  on  the  time-history  of  the  exit  Mach  number  obtained  from  the  difTuser  start-up 
computations  mentioned  above.   Two  major  features  are  noted  in  these  results  (Fig.  4-1). 

The  first  feature  is  the  monotonic  decrease  in  exit  Mach  number  with  time,  as  the  flow  within  the 
difTuser  undergoes  transition  from  an  initial  "cloud  expansion"  mode  to  a  steady  expansion  flow  in  a 
channel  of  increasing  cross-section  area.  A  steady  exit  flow  seems  to  be  established  after  about 
2  (ms)  from  the  start-up  instant.  The  second  feature  is  related  to  the  heat  released  in  the  flow  by 
hydrogen  recombination.    As  the  recombination  rate  is  increased,  so  does  the  time-asymptotic  value 
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of  the  exit  recombination  fraction.  The  recombination  fractions  in  the  three  cases  (Fig.  4-1)  are  1%, 
7%  and  46%;  the  corresponding  exit  Mach  numbers  are  4.0,  3.69  and  2.81.  This  trend  is  in 
qualitative  agreement  with  the  well  known  phenomenon  of  decrease  in  Mach  number  as  a  result  of 
heat  addition  to  a  steady  supersonic  flow  in  a  channel  of  uniform  cross-section  [11]. 

Estimates  of  recombination  rates  are  notoriously  uncertain.  How  can  a  reasonable  uncertainty 
factor  be  established?  We  do  not  know  of  an  estimate  of  this  factor  for  the  flow  of  HF  DF  laser 
products.  However,  in  another  case  of  flow  involving  hydrogen  recombination  -  that  of  hydrazine 
rocket  motors  -  a  recent  study  [12]  recommends  a  hydrogen  recombination  rate  of  2.8  x  1017/T  (in 
c.g.s.  K  units)  and  an  uncertainty  factor  of  30.  Assuming  the  relatively  low  temperature  of  T=300 
(K)  and  multiplying  by  30,  we  get  a  recombination  rate  of  3  x  1016,  versus  7  x  1016  in  our  tenfold 
case.  This  demonstrates  that  the  tenfold  case  is  already  reasonably  high;  just  the  same,  we  also  retain 
the  hundredfold  case  in  the  upcoming  discussion. 

The  exit  Mach  number  time-histories  (Fig.  4-1)  demonstrate  that  the  nominal  difluser  flow  is 
nearly  frozen  (exit  recombination  fraction  1%  and  Mach  number  4);  even  in  the  tenfold  case  the  flow 
is  not  dominated  by  hydrogen  recombination  (exit  recombination  fraction  7%  and  Mach  number 
3.69).  This  conclusion  is  in  agreement  with  findings  of  the  TRW  HF/DF  laser  study  [5]  .  It  takes  the 
unrealistically  high  hundredfold  increase  in  hydrogen  recombination  rate  to  produce  a  decisive  change 
in  the  steady  difluser  flow  (exit  recombination  fraction  46%  and  Mach  number  2.81). 

One  role  of  the  difluser  is  to  expand  the  laser  exhaust  to  a  sufficiently  high  exit  Mach  number,  so 
that  the  thermally  backscattered  flux  of  corrosive  molecules  (HF,  DF)  is  negligibly  small.  Let  us 
consider  the  nominal  case  (M1  =  4)  and  denote  by  "reference  flux"  the  combined  HF+-DF 
backscattered  flux  arriving  at  a  point  located  at  0.1  (m)  from  the  nozzle  lip  (Fig.  1-1).  Using  the  code 
RINGBD  based  on  the  breakdown  surface  model  [2]  ,  the  nominal  reference  flux  was  evaluated  as 
3.1  x  io5  (molecules  m"2  sec"1).  This  flux  level  corresponds  to  a  surface  deposition  rate  (assuming  the 
sticking  coefficient  equals  1)  of  about  10"10  molecular  monolayers  per  hour,  which  is  utterly  negligible 
(a  reasonable  estimate  of  total  operating  time  would  not  exceed  several  hours). 

Observing  that  a  steady  exit  flow  has  been  established  in  all  three  cases  by  about  t  =  2  (ms)  (see 
Fig.  4-1),  we  suggest  that  an  approximate  estimate  of  the  backscattered  flux  would  be  obtained  by 
assuming  a  steady  flow  at  the  difluser  exit  and  the  exhaust  plume.  The  RINGBD  computations  of 
the  increase  hydrogen  recombination  rate  cases  yielded  the  following  results.  In  the  tenfold  case 
(M.  =  3.69)  the  flux  was  102  times  larger  than  the  nominal  reference  flux.    In  the  hundredfold  case 
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(M j  =  2.81)  the  flux  increase  was  by  a  factor  of  107.    Even  in  the  unrealistically  high  hundredfold  case, 
that  flux  level  corresponds  to  about  10"3  monolayers  per  hour,  which  is  still  negligible. 

What  other  uncertainties  could  affect  the  foregoing  conclusion?  First  we  consider  the  effect  of 
variations  in  stagnation  temperature  and  density  relative  to  the  nominal  case,  then  we  take  up  the 
matter  of  transient  plume  flow. 

Exit  stagnation  temperature  was  observed  to  vary  in  the  HFL  diffuser  start-up  computations  by  a 
factor  of  no  more  than  about  1.3.  Since  all  velocities  in  the  breakdown  surface  model  are  normalized 
by  speed  of  sound  [2]  ,  the  flux  is  proportional  to  the  square  root  of  the  stagnation  temperature. 
Also,  it  has  been  shown  [2]  that  a  change  in  stagnation  density  can  produce  an  inversely  proportional 
change  in  flux.  Stagnation  density  was  reduced  by  a  factor  of  about  1.3  in  the  tenfold  case  and  about 
2.8  in  the  hundredfold  case.  In  either  case  the  potential  effect  on  flux  is  small  relative  to  changes 
brought  upon  by  reduced  exit  Mach  number,  so  we  attribute  the  drastic  change  in  flux  following  an 
increased  hydrogen  recombination  rate  primarily  to  the  reduced  exit  Mach  number. 

We  now  address  the  effect  of  transient  plume  flow.  Since  the  exit  Mach  number  is  monotonically 
decreasing  with  time  (Fig.  4-1),  and  since  other  flow  variables  do  not  appreciably  affect  the  flux  from 
a  steady  plume,  we  suggest  that  an  upper  bound  on  the  effect  of  transient  plume  flow  can  be 
established  by  linking  it  to  a  "transient  turning  angle",  larger  than  that  corresponding  to  a  steady  flow 
with  the  same  exit  Mach  number.  Obviously,  if  a  Prandtl-Meyer  flow  pattern  is  derived  from  a  flow 
which  exits  the  nozzle  at  an  angle  lower  than  90°,  its  limiting  turning  angle  would  bring  it  closer  to 
the  spacecraft  surface,  and  with  a  lower  outward  radial  velocity  to  overcome,  more  molecules  would 
be  thermally  backscattered.  Thus,  we  regard  the  transient  flow  around  the  nozzle  lip  as  represented 
by  an  equal  exit  Mach  number  steady  flow,  with  some  initial  exit  turning  angle. 

The  justification  for  this  model  is  the  following  "upper  bound"  estimate  for  the  transient  turning 
angle  of  an  emerging  exhaust  plume.  The  exit  velocity  vector  U(  is  visualized  as  rotating  by  a  gradual 
turning  of  the  plume/ vacuum  interface  until  it  reaches  the  normal  velocity  (2/(y-l))Cj  (Fig.  4-2). 
which  corresponds  to  a  plane-wave  expansion  into  vacuum.  The  corresponding  turning  angle  is 
a  =  (2/(y-l))C]/U1  =  (2/(y-1))/Mj.  It  is  an  overestimate  of  the  turning  angle  since  in  a  steady  corner 
expansion  (and  also  in  a  nonsteady  one)  the  magnitude  of  the  velocity  vector  increases  as  it  rotates, 
in  order  to  conserve  stagnation  enthalpy  (or  even  increase  it  in  nonsteady  flow  about  an  expansive 
corner). 
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It  is  of  interest  to  note  that  a  as  given  above  is  identical  with  the  first  term  of  a  power  series 
expansion  in  1/M,  of  the  standard  Prandtl-. Meyer  expression  for  the  limiting  turning  angle  [13]  .  Let 
us  denote  by  Aa  the  increment  of  the  unsteady  turning  angle  relative  to  the  corresponding  steady  one 
(Prandtl-Meyer).    This  (positive)  increment  is  given  by  : 

Aa  -  (2/(Y-l))/M,   -   {  r,/2arctan[r1/2(M12-l)"1/2]   -  arctan[(M,2  -  l)"'/2]|  } 

(4.1-1) 

r  =  (y+i)/(Y-i) 


We  now  argue  that  the  flux  ratio  between  a  steady  flow  and  the  same  flow  with  exit  velocity 
rotated  by  Aa  is  an  upper-bound  estimate  to  the  effect  of  transient  plume  flow  with  the  same  exit 
variables.  We  note  that  by  pre-rotating  the  steady  corner  How,  we  make  the  limiting  streamline  in  the 
steady  case,  coincide  with  the  nonsteady  plume/vacuum  interface  deemed  to  have  rotated  through  the 
angle  a. 

The  stage  is  now  set  to  estimate  the  transient  flux  increase  factors  in  the  three  hydrogen 
recombination  rate  cases  (exit  Mach  numbers  4,  3.69  and  2.81).  Using  Eq.(4-1)  we  get  the  turning 
angle  increments  of  Aa=  4.1°,  5.1°  and  10.6°  respectively.  Re-computing  the  reference  flux  (code 
RIXGBD)  in  these  cases  with  initial  exit  angle  of  90°- Aa,  we  get  the  flux  increase  factors  of  2.7.  3.5 
and  10  respectively.  Even  in  the  worst  case  (hundredfold),  the  flux  would  increase  from  10"3  to  10"" 
monolayers  per  hour,  which  is  still  negligibly  small. 

We  also  notice  that  the  flux  increase  factor  resulting  from  the  decrease  in  exit  Mach  number  due 
to  higher  hydrogen  recombination  rate,  is  much  larger  than  the  factor  representing  the  effect  of 
transient  plume  flow  in  the  same  case.  Thus,  in  the  tenfold  case  these  factors  are  10"  versus  3.5,  and 
in  the  hundredfold  case  the  ratio  is  even  higher  :  107  versus  10.  We  conclude  that  the  transient  effect 
is  quantitatively  secondary  to  the  effect  of  exit  Mach  number,  which  is  thus  established  as  the  major 
parameter  for  designing  an  exhaust  with  a  negligible  level  of  thermally  backscattered  HF+DF  flux. 
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4.2      Kinetic   Backscattering 

The  role  of  the  third  body  in  the  hydrogen  recombination  reaction  (3.2-1)  is  to  absorb  the  energy 
released  by  the  recombination  process  while  maintaining  the  combined  pre-collision  momentum  of 
the  participating  molecules.  Assume  the  momentum  added  to  [M]  is  q,  then  the  momentum  added  to 
H,  is  -q.  Denote  by  m{  the  mass  of  H2,  by  ra,  the  mass  of  HF  and  by  e  the  energy  released  due  to 
the  recombination.    Then  q  is  determined  from  the  following  energy  conservation  relation  : 

q2/2m,   +  q2/2m2  =  e 

(4.2-1) 
m,  =  2  m2  =  20  e  =  4.36  x  108  (Joules  kmole) 

U  =  q/m,  =  2000   (msec) 

The  values  used  in  the  computation  above  are  per  kmole,  but  the  result  is  the  same  as  using  values 
pertaining  to  single  molecules.  Since  the  velocity  obtained  upon  expansion  to  zero-pressure  in  the 
typical  case  is  about  3000  (m/sec),  our  results  imply  that  a  turning  angle  in  excess  of  arccos(2;3)  =  48° 
is  needed  in  order  to  enable  some  kinetically  backscattered  molecules  to  reach  the  spacecraft. 

In  the  nominal  case  the  turning  angle  is  49°  and  in  the  tenfold  case  it  is  52°.  Also,  some 
kinetically  backscattered  molecules  can  originate  from  within  the  rarefaction  fan,  since  the  exit 
velocity  in  the  typical  case  is  only  about  1500  (m  sec"1)  (versus  2000  (m  sec"1)  velocity  increment  to 
kinetically  scattered  HF  molecules).  Consequently,  this  effect  may  contribute  to  spacecraft 
contamination,  and  its  magnitude  should  be  estimated.  It  will  become  progressively  more  significant 
at  lower  exit  Mach  number,  and  thus  may  be  an  additional  factor  in  determining  a  minimal  value  of 
this  flow  variable. 

While  we  do  not  presently  attempt  at  an  exact  integration  of  the  flux  of  kinetically  backscattered 
molecules  arriving  at  each  surface  point,  we  propose  the  following  overestimate  of  its  magnitude. 
Consider  the  outgoing  flux  of  kinetically  scattered  HF  molecules  from  a  half-space  of  quiescent 
uniform  gas  having  the  thermodynamic  properties  of  the  nozzle  exit  in  the  typical  case.  Consider  the 
hydrogen  recombination  rate  equation  (3.2-1),  and  substitute  [X]  =  fxn  Ay.  This  equation  then  reads 
as  the  volume  rate  of  hydrogen  recombinations  (denoted  as  g  below).  In  a  slab  of  quiescent  gas,  half 
the  kinetically  scattered  molecules  have  an  outward  pointing  velocity  vector  in  one  direction.  The 
probability  of  collisionless  passage  out  of  the  slab  is  exp(  —  x/X),  so  that  the  outgoing  flux  F  is  given 
by: 
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F  -  gc  ^/2 

(4.2-2) 
g    =  kft  n3  A  "2  f  2  (ful7  +  f„c) 

"c  0  v       H     v  HF  DF7 

Using  typical  case  exit  conditions  n,  =  2.81  x  1022  (m"3)  and  mean  free  path  X,  =  1.28  x  10"4  (m), 
we  get  F=5.14*  10  (m""  sec"1),  or  about  12  monolayers  per  hour.  Since  only  a  small  fraction  of 
this  flux  actually  arrives  at  spacecraft  surface  points,  the  arriving  flux  is  in  all  likelihood  well  below  1 
monolayer  per  hour.  If  this  flux  level  is  deemed  significant,  an  accurate  integration  of  kinetically 
scattered  flux  emanating  from  the  exhaust  plume  should  be  performed.  The  details  of  such  scheme 
are  analogous  to  the  first-collision  ambient  scattering  model  [4]  ,  since  both  effects  result  in  a  source- 
like distribution  throughout  the  plume,  and  both  involve  a  steric  factor  (one  half  in  the  quiescent  gas 
above)  related  to  the  vector  addition  of  flow  velocity  and  velocity  imparted  through  ambient  collision 
or  assistance  in  hydrogen  recombination.  Therefore,  the  two  effects  can  conveniently  be  unified 
under  a  common  framework. 


4.3      Concluding   Remarks 

From  the  foregoing  analysis  and  discussion  we  conclude  that  a  design  for  low  (negligible)  level  of 
thermally  backscattered  flux  of  HF+DF  can  be  accomplished  by  assuming  steady  frozen  exhaust 
flow.  The  effects  of  transient  difiuser,  plume  flow  and  hydrogen  recombination  may  typically  increase 
that  flux  by  a  factor  of  102.  Given  the  sensitivity  of  flux  to  exit  Mach  number,  an  adequate  margin 
can  readily  be  achieved  by  raising  the  exit  Mach  number.  In  the  typical  case  presented  here,  an  exit 
Mach  number  of  4  has  been  shown  to  be  adequate. 

It  should  be  noted  that  transient  effects  originating  from  points  upstream  of  the  diffuser  entrance 
were  not  considered  here,  as  they  depend  on  details  of  the  system  construction  and  operating 
sequence. 

The  kinetic  scattering  effect  has  been  pointed  out  as  an  additional  potential  source  of 
contaminating  backflow.  It  is  the  result  of  hydrogen  recombination  collisions  assisted  by  an  IIF  or 
DF  molecule  as  third  body.  The  HF  DF  molecule  receives  part  of  the  kinetic  energy  release  by  the 
hvdrosen  recombination. 
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APPENDIX   A.        The   HFL   Code 

A.  1     Features  Specific   to   Diffuser   Flow   Model 

The  purpose  of  this  appendix  is  to  provide  the  listing  of  the  HFL  code  used  in  the  diffuser  How 
computations.  An  almost  identical  code  was  recently  described  in  some  detail  [7]  ,  the  major 
modification  for  HFL  being  the  addition  of  hydrogen  recombination. 

The  highlights  of  the  HFL  version  to  the  GRP  code  [7]  are  as  follows  (numbers  in  parenthesis 
refer  to  statement  numbers  which  are  prefixed  by  HFL  in  the  listing). 

(a)  Data  (NETUNM)  : 

The  data  which  relates  specifically  to  the  diffuser  flow  with  hydrogen  recombination  is  given  in 
statements  (194  -  239).  Some  parameters  are  stored  in  /DIFFUS/  (129)  for  use  in  BEGIN  and 
SAFAE. 

(b)  Initial  Conditions  (BEGIN)  : 

The  initial  conditions  are  ideally  an  empty  diffuser.  As  an  approximation  we  set  small  pressure 
and  density  (270  -  271)  and  a  positive  velocity  (277)  as  initial  values.  The  velocity  reduces  the 
intensity  of  the  shock  driven  into  the  dilute  initial  gas  by  the  abrupt  start-up  of  inflow  at  the 
diffuser  entrance  (1  =  2).  Note  that  the  initial  value  of  the  total  energy  E(I)  is  augmented  by  the 
heat  of  formation  of  hydrogen  (306).  The  heat  term  is  subsequently  subtracted  from  the  total 
energy  in  order  to  compute  the  pressure  (see  (d)  below). 

The  grid  is  defined  as  follows.  The  number  of  cells  (L-2)  is  L-2=50  (24).  The  diffuser  extends 
from  X0  =  0.625  to  XI  =  3.125  (283  -  284),  so  that  DX  =  0.05.  Since  the  external  radius  of  the 
spacecraft  is  just  2.5,  this  grid  consists  of  40  cells  within  the  diffuser  and  10  cells  in  an  extended 
segment.  The  purpose  of  this  extended  segment  is  to  minimize  any  propagation  of  error  due  to 
imperfection  of  the  outflow  boundary  condition  (see  (c)  below).  Since  the  ambient  pressure  is 
zero,  this  additional  segment  does  not  affect  the  flow  at  X  £  2.5.  Consequently,  the  results  for 
the  diffuser  exit  flow  were  read  from  cell  I  =  39  whose  midpoint  is  X  =  2.5. 

(c)  Boundary  Conditions  (SAFAE)  : 

The  inflow  is  set  by  assigning  the  diffuser  entrance  conditions  to  cell  1=1.  The  outflow  is  set 
by  extrapolating  (with  zero  gradient)  the  flow  from  the  inner  boundary  cell  (I  =  L-1)  to  the 
boundary  cell  I  =  L. 
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(d)      Conservation  Laws   (CYCEUL)  : 

Here  all  four  conservation  laws  are  integrated  in  time.  There  is  an  added  equation  of  a  "time- 
split"  scheme  for  solving  the  conservation  of  H  species  equation.  It  is  first  solved  without  the 
recombination  rate  term  (459),  then  the  change  in  Z(I)  due  to  recombination  during  that  time 
step  is  added  as  DZZ  (501  -  502). 

The  energy  equation  includes  the  contribution  of  the  formation  energy  indirectly.  This  equation 
has  the  same  format  as  the  adiabatic  one  (491),  but  when  QDET.GT.O,  the  internal  energy 
(505)  is  affected  by  changes  in  Z(I)  and  as  a  result  the  pressure  P(I)  (516)  is  also  affected  by- 
progress  in  the  recombination  reaction. 
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A. 2     Code   Listing 

C$0PTI0NS  LIST' -  HFL0001 

IMPLICIT  REAL*8CA-H,0-Z,$)  HFL0002 

C   PROGRAM  DETO  HFL0003 

C   HFL  --   HF/DF  LASER  DIFUSER  TRANSIENT  FLOW.  HFL0004 

COMMON  B(52,26)  HFL0005 

1       ,ENDB  HFL0006 

COMMON  /AB/A(50)  HFL0007 

EQUIVALENCE  ( L , AC  1 ) ) , ( LL , A( 2) ) , ( T, AC  3) ) , ( DT, A( 4) ) , CTMAX, A( 5 ) ) ,     HFL 0  0  08 

1  CTMUD,AC6)), CDTMUD,AC7)), CJ0B,A(8)),CNERI,AC9)),  HFL0009 

2  CJJJ,AC10)),  CKEYM0N,AC11)),CNCYC,AC12))  HFL0010 
EQUIVALENCE  CCOLELA, AC  13) )  HFLOOll 
EQUIVALENCE  C LAGEUL , AC  14 ) )  HFL0012 
EQUIVALENCE  CUGAL,AC15))  HFL0013 
EQUIVALENCE  C KEYEK, AC  16 ) )  HFL0014 
EQUIVALENCE  C NCYCPR, AC  17 ) )  HFL0015 
EQUIVALENCE  C STAB, AC  18 ) ) , C DTBA, AC  19 ) ) , C DTKOD, AC  20) ) , C KDT, AC  21) )  HFL 00 16 
COMMON  /M0NIT/NC14C4),CASEAVC4),NF16C6),  HFL0017 

1                NM0NUC4),NM0NPC4),NM0NGC4),NM0NR0C4),NM0NZC4)  HFL0018 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXHFLOOig 

C      DO  20  N=l,30  HFL0020 

C20    NC14CN)=0  HFL0021 

NMAT=26  HFL0022 

C      L=CL0CFCENDB)-L0CFCBC1,1)))/NMAT  HFL0023 

L=52  HFL0024 

LL=L-1  HFL0025 

NN=NMAT*L  HFL0026 

DO  1  1=1, L  HFL0027 

DO  1  11=1, NMAT  HFL0028 

1     BCI,II)=0.  HFL0029 

CALL  MAIN0CL,BC1,  1),BC1,  2),BC1,  3),BC1,  4),BC1,  5),  HFL0030 

1  BC1,  6),BC1,  7),BC1,  8),BC1,  9),BC1,10),  HFL0031 

2  BC1,11),BC1,12),BC1,13),BC1,14),BC1,15),  HFL0032 

3  BC1,16),BC1,17),BC1,18),BC1,19),BC1,20),  HFL0033 

4  BC1,21),BC1,22),BC1,23),BC1,24),BC1,25),  HFL0034 

5  BC1,26))  HFL0035 
STOP  HFL0036 

< EMT2 HFL0037 

SUBROUTINE  MAINO                             "  HFL0038 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0039 

2  US,PS,UIDOT,PIDOT,  HFL0040 
X             FIMZ,ZMDOT,  HFL0041 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0042 
IMPLICIT  REAL*8CA-H,0-Z,$)  HFL0043 
DIMENSION  XCL),UCL),PCL),ROCL),GCL),ECL),DUCL),DPCL),DROCL),  HFL 0  044 

1  DGCL),DXSICL),MINCL),  HFL0045 

2  USCL),PSCL),UIDOTCL),PIDOTCL)  HFL0046 

3  ,TENACL),FIROCL),FIMCL),FIECL)  HFL0047 

4  ,GIPCL),VOLCL),ZCL),DZCL)  HFL0048 

5  ,FIMZCL),ZMDOTCL)  HFL0049 
COMMON  /AB/AC50)  HFL0050 
EQUIVALENCE  C LL , AC2) ) , CT, AC  3) ) , C DT, AC  4) ) , CTMAX, AC  5) ) ,  HFL0051 

1  CTMUD,AC6)),CDTMUD,AC7)),CJ0B,AC8)),CNERI,AC9)),  HFL0052 

2  CJJJ,AC10)),CKEYMON,AC11)),CNCYC,AC12))  HFL0053 
EQUIVALENCE  C LAGEUL , AC  14) )  HFL0054 
EQUIVALENCE  C NCYCPR, AC  17  ) )  HFL0055 
EQUIVALENCE  C STAB, AC  18 ) ) , C DTBA, AC  19 ) ) , C DTKOD, AC  20 ) ) , CKDT, AC 21 ) )  HFL 0  056 
COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  HFL0057 

CXXXXXXXXXXXX3^(XXXXXXXX3(XXXXXXXXXXXXXXXXXKKKXXXXXXXXXKX*XXXXXKXXX*X)(XXXXHFL0  058 

T=0.  HFL0059 

NCYC=0  HFL0060 

JJJ=0  HFL0061 

CALL  NETUNM  HFL0062 

DELT=DT  HFL0063 

CALL  BEGIN  HFL0064 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0065 

2  US,PS,UIDOT,PIDOT,  HFL0066 
x             FIMZ,ZMDOT,  HFL0067 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0068 
CALL  SAFAE  HFL0069 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0070 

2  US,PS,UIDOT,PIDOT,  HFL0071 
x             FIMZ,ZMDOT,  HFL0072 
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LE:  HFL 


SCRIPT 


Al 


3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0073 

1  NCYC  =  NCYC+1  HFL0074 
C  TIME  STEP  CONTROL.  HFL0075 

DT=DTBA"  HFL0076 

IFCDT.GT.l .1D0XDTKOD.AND.DTKOD.NE.0.  )  DT=1 . 1DO*DTKOD  HFL0077 

IF(NCYC.EQ.2)  DT=DT/10.  HFL0078 

IF  (NCYC.EQ.l)  DT=0.  HFL0079 

IFCDT.EQ.O.)  GO  TO  11  HFL0080 

NHAD=((TMUD-T)/DT-1.D-10)  HFL 0  081 

IF(NHAD.GE.IO)  GO  TO  11  HFL0082 

DT=(TMUD-T)/DFL0AT(NHAD+1)  HFL 0  08  3 

11    CONTINUE  HFL0034 

T=T+DT  HFL0085 

IF((NCYC/NCYCPR)*NCYCPR.NE.NCYC.AND.NCYC.GT.NCYCPR)  GO  TO  33  HFL0086 

PRINT  10,  NCYC,T,DT,KDT  HFL0087 
10    FORMAT (IX, ,NCYC=',I^,3X, 'T= ' , Dll . 4, 3X, ' DT= ' , Dll . 4, 3X, fKDT=M4)    HFL 0088 

33    CONTINUE  HFL0089 

DTBA=DTMUD  HFL0090 

KDT=0  HFL0091 

NERI=1  HFL0092 

IF  (DABS(T-TMUD).LT.l.D-8)  NERI=0  HFL0093 

CALL  CYCEUL  HFL0094 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0095 

2  US,PS,UIDOT,PIDOT,  HFL0096 

*  FIMZ,ZMDOT,  HFL0097 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0098 
CALL  SAFAE  HFL0099 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0100 

2  US,PS,UIDOT,PIDOT,  HFL0101 

*  FIMZ,ZMDOT,  HFL0102 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0103 
IF  (NERI.NE.O)  GO  TO  2  HFL0104 
CALL  PRINT  HFL0105 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0106 

2  US,PS,UIDOT,PIDOT,  HFL0107 
X             FIMZ,ZMDOT,  HFL0108 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0109 
IF  (DABS(T-TMUD) .LT.l.D-8)  TMUD=TMUD+DTMUD  HFL0110 

2  CONTINUE  HFL0111 
DTKOD=DT  HFL0112 
IF  (T.LT.TMAX-l.D-8)  GO  TO  1  HFL0113 
RETURN  HFL011A 
END HFL  0115 

SUBROUTINE  NETUNM RTTOTT^ 

IMPLICIT  REAL*8(A-H,0-Z,$)  HFL0117 

COMMON  /AB/A(50)  HFL0118 

EQUIVALENCE  (L,A(D)  HFL0119 

EQUIVALENCE  ( LL , A(2) ) , (T, A( 3) ) , ( DT, A(4) ) , (TMAX, A(5) ) ,  HFL0120 

1  (TMUD,A(6)),(DTMUD,A(7)),(J0B,A(8)),(NERI,A(9)),  HFL0121 

2  (JJJ,A(10)),(KEYMON,A(11)),(NCYC,A(12))  HFL0122 
EQUIVALENCE  (COLELA, A( 13) )  HFL0123 
EQUIVALENCE  ( LAGEUL , A( 14) )  HFL0124 
EQUIVALENCE  ( KEYEK, A( 16  ) )  HFL0125 
EQUIVALENCE  ( NCYCPR, A( 17 ) )  HFL0126 
EQUIVALENCE  ( STAB, A( 18 ) ) , ( DTBA, A( 19) ) , ( DTKOD, A(20) ) , CKDT,  A(  21 ) )  HFL 01 27 
COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET  HFL 01 28 
C0MM0N/DIFFUS/U2,P2,R02,ARW  HFL 01 29 
COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX  HFL 01 30 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA  HFL0131 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 01 32 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20  HFL0133 

REALX8  NG  HFL0134 

REAL*8  MU2  HFL0135 

NAMELIST  /IN/LIN, GAMA, DT, TMUD, DTMUD, TMAX,  HFL0136 

1  GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX,  HFL0137 

2  SMIN,SMAX,IVERSA,KEYMON, COLELA, STAB  HFL0138 

3  ,LAGEUL, KEYEK  HFL0139 

4  ,QDET,TC,RATE  HFL0140 

CXXXXXXXXXXXXXXXXXXXXXXXXXX^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXHFLOl^l 

LIN=L  HFL0142 

LAGEUL=2  HFL0143 

NCYCPR=20  HFL0144 
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HFL       SCRIPT    Al 

KEYEK=1  HFL0145 

TMUD=0.  HFL0146 

DTMUD=0.2D0  HFL0147 

TMAX=2.0D0  HFL0148 

DT=2.5D-3  HFL0149 

GAMA=1.54D0  HFL0150 

KEYM0N=1  HFL0151 

STAB=0.6D0  HFL0152 

C      RATE=1.5D0  HFL0153 

G0DELX=16.  HFL0154 

GODELY=20.  HFL0155 

IVERSA=100  HFL0156 

UMIN=0.  HFL0157 

UMAX=  3. DO  HFL0158 

PMIN=0.  HFL0159 

PMAX=l.D-3  HFL0160 

ROMIN=0.  HFL0161 

R0MAX=l.D-3  HFL0162 

SMIN=0.  HFL0163 

SMAX=PMAX/ROMAXXXGAMA  HFL0164 

COLELA=0.  HFL0165 

C      READ  IN  HFL0166 

C  HFL0167 

PRINT  IN  HFL0168 

C  HFL0169 

GG=2.*GAMA/(GAMA-1. )  HFL0170 

NG=GG  HFL0171 

10    CONTINUE  HFL0172 

MU2=(GAMA-1.)/(GAMA+1.)  HFL0173 

G1=GAMA-1.  HFL0174 

G2=1.-MU2  HFL0175 

G3=2./(3.XGAMA-1.)  HFL0176 

G4=(GAMA+l.)/2.  HFL0177 

G5=0.5*(3.*GAMA-1. )/(GAMA+l. )  HFL 0178 

G6=(GAMA+1. )/C2.xGAMA)  HFL0179 

G7=2./(GAMA-1.)  HFL0180 

G8=(GAMA-1.)/(2.XGAMA)  HFL0181 

G9=(GAMA+1.)/(4.XGAMA)  HFL0182 

G10=1./GAMA  HFL0183 

Gll=(GAMA+l.)/4.  HFL0184 

G12=GAMA/(GAMA-1. )  HFL0185 

G13=0.5x(GAMA-3. )/(GAMA+l.)  HFL 0186 

G14=0.5X(3.XGAMA-5.)/(GAMA+l.)  HFL 0187 

G17=GAMA+1.  HFL0188 

G0DELX=G0DELX/2.54  HFL0189 

G0DELY=G0DELY/2.54  HFL0190 

C      CALL  NAMPLT(IVERSA)  HFL0191 

C      CALL  LIMITC1000.)  HFL0192 

C      CALL  PLOTCO. ,0.5,-3)  HFL0193 

C  DATA  FOR  DIFFUSER  ENTRY  CONDITIONS.  TEST  III/2,  REF.(M-l).  HFL0194 

C  UNITS  ARE  IN  M,K,MILISEC.  HFL0195 

C   FH2  IS  THE  DIFFUSER  ENTRY  MOLAR  FRACTION  OF  H  ATOMS.  HFL0196 

FH2=0.091D0  HFL0197 

U2=2.3D0  HFL0198 

AMD0T=1.499D-5  HFL0199 

A2=1.D0X7.D0*0.0254XX2  HFL0200 

R02=AMD0T/(A2XU2)  HFL0201 

TEMP0=1400.D0  HFL0202 

W2=7.27D0  HFL0203 

ARW=8.31441D-3/W2  HFL0204 

DELTAQ=435. 783/2. DO  HFL0205 

QDET=FH2XDELTAQ/W2  HFL0206 

AA=U2/DSQRT(GAMAXARWXTEMP0)  HFL 0207 

RATE=1 .4D0X1 .D7XFH2/W2XX2  HFL0208 

RATE=RATExi.D2  HFL0209 

EM2=AA/DSQRT(1 .-AAXX2/G7)  HFL0210 

GOREM=(l.D0+EM2**2/G7)  HFL0211 

TRATIO=2.DOX(GAMA+1.DO)XEM2XX2XGOREM/(1.DO+GAMAXEM2XX2)XX2  HFL 0212 

HH0=GAMAXARWXTEMP0/(GAMA-1 .DO)  HFL0213 

HHl  =  HH0/TRATIO  HFL0214 

CHOKE=QDET/(HH1-HHO)  HFL0215 

TEMP2=TEMP0/G0REM  HFL0216 
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LE=  HFL       SCRIPT    Al 

P2=ARW*R02*TEMP2  HFL0217 

RATIO=2.5D0/0.5D0  HFL0218 

EM0UT=EM2*RATI0xx((GAMA-l  . D0)/2.D0)  HFL0219 

DO  250  1=1,100  HFL0220 

EM0UT1=EM0UT  HFL0221 
GOREM=RATIO*EMOUTX(1.DO+0.5DOX(GAMA-1.DO)*EM2X*2)*X(0.5DO/MU2)/EM2HFL0  222 

G0REM=GOREMX*(2.D0XMU2)  HFL0223 

EMOUT=DSQRT((2.DO/(GAMA-1 . DO ) ) x(GOREM-l .DO))  HFL 0224 

DM=DABS(EM0UT-EM0UT1)  HFL0225 

IFCDM.LT.l.D-9)  GO  TO  251  HFL0226 

250  CONTINUE  HFL0227 
CALL  S0F(251)  HFL0223 

251  CONTINUE  HFL0229 
PRINT  201  HFL0230 

201  F0RMAT(//1X, 'DIFFUSER  ENTRY  DATA')  HFL0231 
PRINT  202,U2,P2,R02,FH2  HFL0232 

202  F0RMAT(/1X,  »U2, P2, R02, FH2= ' , 4D16 . 5)  HFL0233 
PRINT  203, AMD0T,W2, TEMPO, TEMP2,EM2  HFL0234 

203  F0RMATC/1X, ' AMDOT, W2, TEMPO, TEMP2, EM2= ' , 5D16 . 5)  HFL0235 
PRINT  204, RATE, QDET, HHO , HH1 , CHOKE  HFL0236 

204  F0RMAK/1X,  '  RATE,  QDET,  HHO  ,  HH1 ,  CHOKE=  '  ,  5D16  .  5/ )  HFL0237 
PRINT  205, RATIO, EMOUT  HFL0238 

205  F0RMATC1X, 'EXIT  AREA  RATIO  AND  MACH  NUMBER   RATIO, EMOUT= ' , 2D16 . 5)  HFL0239 
RETURN  HFL0240 

END    HFLQ241 

SUBROUTINE  BEGIN HFL0242 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0243 

2  US,PS,UIDOT,PIDOT,  HFL0244 
X             FIMZ,ZMDOT,  HFL0245 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0246 
IMPLICIT  REALX8(A-H,0-Z,$)  HFL0247 
DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  HFL 0248 

1  DG(L),DXSI(L),MIN(L),  HFL0249 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  HFL0250 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  HFL0251 

4  ,GIP(L),VOL(L),Z(L),DZ(L)  HFL0252 

5  ,FIMZ(L),ZMDOT(L)  HFL0253 
COMMON  /AB/AC50)  HFL0254 
COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 0255 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20  HFL0256 

REALX8  NG  HFL0257 

REALX8  MU2  HFL0258 

EQUIVALENCE  (LL,A(2))  HFL0259 

EQUIVALENCE  ( LAGEUL , A( 14)  )  HFL0260 

EQUIVALENCE  (UGAL,A(15))  HFL0261 
EQUIVALENCE  ( STAB, A( 18 ) ) , ( DTBA, A( 19 ) ) , ( DTKOD, A (20 ) ) , (KDT, A( 21 ) )  HFL 0  26 2 
COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET  HFL 026 3 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX  HFL 026  4 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA  HFL0265 

CXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXHFL0266 

DTBA=0.  HFL0267 

DTKOD=0.  HFL0268 

KDT=0  HFL0269 

P0=l.D-7  HFL0270 

RHO0=l.D-7  HFL0271 

Pl=0.1  HFL0272 

Pl=4.  HFL0273 

RH01=0.125  HFL0274 

RH01=4.  HFL0275 

UGAL=0.  HFL0276 

U0=1.08D0  HFL0277 

U1=0.  HFL0278 

UO=UO-UGAL  HFL0279 

U1=U1-UGAL  HFL0280 

UMIN=UMIN-UGAL  HFL0281 

UMAX=UMAX-UGAL  HFL0282 

X0=0.625D0  HFL0283 

X1=3.125D0  HFL0284 

XMIN=XO  HFL0285 

XMAX=X1  HFL0286 

DX=(Xl-X0)/(L-2.)  HFL0287 

DO  1  1=2, L  HFL0288 
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X(I)=X0+(I-2. )*DX  HFL0289 

1  CONTINUE  HFL0290 
X(L)=X1  HFL0291 
XD=1.D10  HFL0292 
DO  2  1=2, LL  HFL0293 
U(I)=UO/X(I)  HFL0294 
P(I)=PO  HFL0295 
RO(I)=RHOO  HFL0296 
Z(I)=0.  HFL0297 
IF  (XCD-XD.LT.-1.D-8)  GO  TO  21  HFL0298 
U(I)=U1  HFL0299 
P(I)=P1  HFL0300 
RO(I)=RH01  HFL0301 
Z(I)=l.DO  HFL0302 

21    CONTINUE  HFL0303 

GO  TO  (31,32),  LAGEUL  HFL0304 

31  CONTINUE  HFL0305 
E(I)=P(I)/((GAMA-1. )XRO(I))+0.5XU(I)X*2+Z(I)XQDET  HFL 0  306 
GO  TO  30  HFL0307 

32  CONTINUE  HFL0308 

E(I)=P(I)/(GAMA-1.)+0.5*RO(I)XU(I)*X2+Z(I)XRO(I)*QDET  HFL 0309 

30    CONTINUE  HFL0310 

G(I)=DSQRT(GAMA*P(I)*RO(I))  HFL 0311 

2  CONTINUE  HFL0312 
DO  3  1=2, LL  HFL0313 
DXSI(I)=(X(I+l)-X(I))xRO(I)  HFL 0314 
TENA(I)=RO(I)*U(I)  HFL0315 
VOL(I)=(X(I+l)-X(I))x(X(I+l)+X(I))/2.D0  HFL 0316 

3  CONTINUE  HFL0317 
RETURN  HFL0318 

END HFLQ319 

DOUBLE  PRECISION  FUNCTION  RATIO(X)  HFL0320 

IMPLICIT  REAL*8(A-H,0-Z,$)  HFL0321 

IF(X.LE.1.D-8)RATIO=0.  HFL0322 

IF(X.LE.1.D-8)RETURN  HFL0323 

RATIO=l.D0/X  HFL0324 

RETURN  HFL0325 

END HFL0326 

DOUBLE  PRECISION  FUNCTION  CR05S(X) H\-LU62/ 

IMPLICIT  REALX8(A-H,0-Z,$)  HFL0328 

CROSS=X  HFL0329 

RETURN  HFL0330 

END HFL0331 

SUBkUUIlNb  CYCEUT HFUTTS2 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL0333 

2  US,PS,UIDOT,PIDOT,  HFL0334 
x             FIMZ,ZMDOT,  HFL0335 

3  TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL0336 
IMPLICIT  REALX8(A-H,0-Z,$)  HFL0337 
DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  HFL 0338 

1  DG(L),DXSI(L),MIN(L),  HFL0339 

2  US(L),PS(L),UIDOT(L),PIDOT(L)  HFL0340 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  HFL0341 

4  ,GIP(L),VOL(L),Z(L),DZ(L)  HFL0342 

5  ,FIMZ(L),ZMDOT(L)  HFL0343 
COMMON  /AB/AC50)  HFL0344 
EQUIVALENCE  ( LL , A( 2) ) , (T, A( 3) ) , ( DT, A(4) ) , (COLELA, A( 13) )  HFL0345 
EQUIVALENCE  ( KEYEK, A( 16 ) )  HFL0346 
EQUIVALENCE  (STAB, A( 18) ) , ( DTBA, A( 19 ) ) , ( DTKOD, A (20 ) ) , ( KDT, A(21 ) )  HFL 0  347 
COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 0  348 

1           ,G12,G13,G14,G15,G16,G17,G18,G19,G20  HFL0349 

REALX8  NG  HFL0350 

REALX8  MU2  HFL0351 

COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  HFL0352 

COMMON  /AZOV/ISAFA,NORIMN,USAF,PSAF,ROSAF,GSAF,ESAF,DPSAF  HFL 0  353 

1             ,DXSIL,DXSIR  HFL0354 

LOGICAL  NORIMN  HFL0355 

COMMON  /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  HFL0356 

1  RSTARL,RSTARR,GSTARL,GSTARR,WL,WR,  HFL0357 

2  CL,CR,CSTARL,CSTARR,UW(6)  HFL0358 

3  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  HFL0359 
LOGICAL  HELEML,HELEMR  HFL0360 
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COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR  HFL0361 

2  ,ASTARL,ASTARR  HFL0362 

3  ,RAT,SH  HFL0363 
COMMON  /GRADS/DUDXIL,DPDXIL,DGDXIL,DRDXIL,  HFL0364 

1              DUDXIR,DPDXIR,DGDXIR,DRDXIR  HFL0565 

COMMON  /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN  HFL0366 

1  ,GIH  HFL0367 

2  ,FIH4,ZMD0TL,ZMD0TR  HFL0368 
COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET  HFL 0  36  9 

C      DATA  K0TZ/7777777777B/  HFL0370 

CXXXXXXXX5<X^XX^X3(^X>(3<3()(XXX3(XXXX)(5(5(XXXXX)(X^XXXXXXXX^X^XX5(XXXXXXXXXXXXX*X^HFL0371 

TN  =  T  HFL0372 

DT2=DT/2.  HFL0373 

DO  3  1=1, L  HFL0374 

C      MINO=MIN(I) .AND.KOTZ  HFL0375 

C      MINC I)=SHIFT(MINO,30)  HFL0376 

3   CONTINUE  HFL0377 

UXN=0.  HFL0378 

PXN=0.  HFL0379 

GXN=0.  HFL0380 

ZN=0.  HFL0381 

ROXN=0.  HFL0382 

DO  1  1=2, L  HFL0383 

IM=I-1  HFL0384 

UXNM=UXN  HFL0385 

PXNM=PXN  HFL0386 

GXNM=GXN  HFL0387 

ROXNM=ROXN  HFL0388 

ZNM=ZN  HFL0389 

UL=U(IM)+0.5XDU(IM)  HFL0390 

PL=P(IM)+0.5*DP(IM)  HFL0391 

ROL=RO(IM)+0.5*DRO(IM)  HFL0392 

GL=G(IM)+0.5*DG(IM)  HFL0393 

CL=GL/ROL  HFL0394 

ZL=Z(IM)+0.5*DZ(IM)  HFL0395 

UR=U(I)-0.5XDU(I)  HFL0396 

PR=P(I)-0.5*DP(I)  HFL0397 

GR=G(I)-0.5*DG(I)  HFL0398 

ROR=RO(I)-0.5*DRO(I)  HFL0399 

CR=GR/ROR  HFL0400 

ZR=Z(I)-0.5XDZ(I)  HFL0401 

CALL  RIEMAN(L,I,MIN)  HFL0402 

DUDXIL  =  DU(IM)/DXSI(IM)  HFL0403 

DPDXIL=DP(IM)/DXSI(IM)  HFL0404 

DGDXIL=DG(IM)/DXSI(IM)  HFL0405 

DRDXIL=DRO(IM)/DXSI(IM)  HFL0406 

DUDXIR  =  DU(I)/DXSI(I)  HFL0407 

DPDXIR=DP(I)/DXSI(I)  HFL0408 

DGDXIR=DG(I)/DXSI(I)  HFL0409 

DRDXIR=DRO(I)/DXSI(I)  HFL0410 

SH  =  CROSS(X(D)  HFL0411 

RAT  =  RATIO(X(D)  HFL0412 

CALL  MAGA(L,I,MIN)  HFL0413 

US(I)=USTAR  HFL0414 

PS(I)=PSTAR  HFL0415 

UIDOT(I)=DUIDT  HFL0416 

PIDOT(I)=DPIDT  HFL0417 

CALL  FLUXE1(L,I,MIN)  HFL0418 

CALL  FLUXE2(L,I,MIN)  HFL0419 

FIR0(I)=FIH1  HFL0420 

FIM  (I)  =  FIH2  *i\™l\ 

FIE  (I)  =  FIH3  H&SJI1 

GIP(I)=GIH  HFL0423 

FIMZ(I)=FIH4  HEh25«= 

DU(IM)=UXN-UXNM  HFL0425 

DP(IM)=PXN-PXNM  HEf-SrSS 

DG(IM)=GXN-GXNM  NFL0427 

DRO(IM)=ROXN-ROXNM  tlEf*  25S5 

1    CONTINUE  SE5-2512 

AMTOT  =  0.  H&SJJ? 

ETOT  =  0.  HFL0431 

EKTOT=0.  HFL0432 
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7352 

7353 
7351 
7350 
2 

201 


EPT0T=0. 

HFL0433 

TENTOT=0. 

HFL0434 

FI1=FIR0(2) 

HFL0435 

FI2=FIM  (2) 

HFL0436 

FI4=FIMZ(2) 

HFL0437 

GI2=GIP(2) 

HFL0433 

SH=CR0SS(X(2)) 

HFL0439 

DO  2  1=2, LL 

HFL0440 

IP=I+1 

HFL0441 

FIM1=FI1 

HFL0442 

FIM2=FI2 

HFL0443 

FIM4=FI4 

HFL0444 

GIM2=GI2 

HFL0445 

SHM=SH 

HFL0446 

GI2=GIP(IP) 

HFL0447 

SH=CR05S(X(IP)) 

HFL0448 

DVOL=VOL(I) 

HFL0449 

FI1=FIR0(IP) 

HFL0450 

FI2=FIM  (IP) 

HFL0451 

FI4=FIMZ(IP) 

HFL0452 

ZKODM=Z(I)XRO(I) 

HFL0453 

DX=X(IP)-X(I) 

HFL0454 

R0(I)=R0(I)-DT/DV0LX(SHXFI1-SHM*FIM1) 

HFL0455 

TENA(I)=TENA(I)-DT/DV0Lx(SHxFI2-SHMxFIM2)-DT/DXX(GI2-GIM2) 

HFL0456 

U(I)=TENA(I)/RO(I) 

HFL0457 

IF(QDET.EQ.O. )  GO  TO  2 

HFL0458 

Z(I)=(ZK0DM-(DT/DV0L)x(SHxFI4-SHMxFIM4))/R0(I) 

HFL0459 

DZZ1=Z(I)-1.D0 

HFL0460 

IF(DZZ1.LT.1.D-6.AND.Z(I) .GE.O.)  GO  TO  7350 

HFL0461 

IFCDZZ1 .LT.O.)  GO  TO  7352 

HFL0462 

Z(I)=1.D0 

HFL0463 

IFCDZZl.LT. 0.2D0)  GO  TO  7350 

HFL0464 

CONTINUE 

HFL0465 

IF(Z(I) .GT.O. )  GO  TO  7353 

HFL0466 

Z(I)=0. 

HFL0467 

IF(DABS(Z(I)).LT.0.2D0)  GO  TO  7350 

HFL0468 

CONTINUE 

HFL0469 

PRINT  7351,I,Z(I),P(I),R0(I),E(I),U(I),ZMD0T(I) 

HFL0470 

FORMAT (/1X, "CYCEUL.  I , Z,  P,  RO, E, U, ZMDOT= • , 15, 6D15 . 5/) 

HFL0471 

CALL  SOFC7350) 

HFL0472 

CONTINUE 

HFL0473 

IF(ZCI) .GT. .1D1)  Z(I)=.1D1 

HFL0474 

CONTINUE 

HFL0475 

IFCCOLELA.EQ.O. )  GO  TO  201 

HFL0476 

CALL  DCOLE(L,X,U  , DU  ,MIN,1) 

HFL0477 

CALL  DC0LE(L,X,R0,DR0,MIN,4) 

HFL0478 

CONTINUE 

HFL0479 

CALL  BD0K1(L,X,U  , DU  ,MIN,1) 

HFL0480 

CALL  BD0K1(L,X,R0,DR0,MIN,4) 

HFL0481 

FI3=FIE  (2) 

HFL0482 

SH=CR0SS(X(2)) 

HFL0483 

DO  6  1=2, LL 

HFL0484 

IP=I+1 

HFL0485 

SHM=SH 

HFL0486 

SH=CROSS(X(IP)) 

HFLQ487 

DVOL=VOL(I) 

HFL0488 

FIM3=FI3 

HFL0489 

FI3=FIE  (IP) 

HFL0490 

DX=X(IP)-X(I) 

HFL0491 

E(I)=E(I)-DT/DV0LX(SHXFI3-SHMXFIM3) 

HFL0492 

HFL0493 

UAV=U(I) 

HFL0494 

HFL0495 

ROAV=RO(I) 

HFL0496 

HFL0497 

EP=E( I )-0 . 5XR0AV*UAVX*2 

HFL0498 

IF  (QDET.EQ.O. )  GO  TO  64 

HFL0499 

TI=P(I) 

HFL0500 

DZZ=-RATEXDTX(R0(I)XZ(I))XX2 

HFL0501 

Z(I)=Z(I)+DZZ 

HFL0502 

IF(Z(I) .LT.O. )  Z(I)=0. 

HFL0503 

IF  (ZCI) .LT.O. 01)  Z(I)=0. 

HFL0504 
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64 
61 
62 


60 


69 
6 


601 


7001 
7101 


EP  =  EP-Z(I)*ROU)xQDET 

CONTINUE 

GO  TO  (61,62),  KEYEK 

CONTINUE 

GO  TO  60 

CONTINUE 

DEK=(R0AV*DU(I)**2+2 

EP=EP-DEK 

GO  TO  60 

CONTINUE 

IFCEP.LE.O. )  GO  TO  7001 

P(I)=G1*EP 

G(I)=DSQRT(GAMAXP(I)*RO(I)) 

DM=RO(I)xDX 

DXSI(I)=DM 

DM=ROCI)xDVOL 

ETOT=ETOT+E(I)XDVOL 

EPTOT=EPTOT+EPxDVOL 

AMTOT=AMTOT+DM 

ETOT=ETOT+E(I)XDX 

EPTOT=EPTOT+EPxDX 

TENTOT=TENTOT+DX*TENA(I) 

UPC=DABS(U(I))+G(I)/RO(I) 

DTI=STABXDX/UPC 

IF(DTI.GT.DTBA)  GO  TO  69 

DTBA=DTI 

KDT  =  I 

CONTINUE 

CONTINUE 

CALL  DCOLE(L 

IFCCOLELA.EQ 

CALL  DCOLECL 

CALL  DCOLECL 

CONTINUE 

CALL  BDOK1CL 

CALL  BDOKKL 

CALL  BDOKKL 

EKTOT=ETOT-EPTOT 

RETURN 

CONTINUE 

PRINT  7101, 

FORMATC//1X, 


xDR0(I)xDU(I)xUAV)/24. 


X,Z,DZ,MIN,5) 
0. )  GO  TO  601 
X,P,DP,MIN,2) 
X,G,DG,MIN,3) 

X,P,DP,MIN,2) 
X,G,DG,MIN,3) 
X,Z,DZ,MIN,5) 


I,ROAV,UAV,DRO(I),DU 
•FROM  CYCEUL .   NEGAT 

1  IX, »ROAV,UAV,DRO(I),DU( 

2  IX, ,E(I),DEK,EP,KEYEK=I 
CALL  SOFC7001) 

RETURN 

FNn 


(I),E( 
IVE  EP 
I)=',4 
,3D18. 


I),DEK,EP,KEYEK 

.   IN  CELL   I=»,I6// 

D18.8// 

8,110//) 


FL0505 
FL0506 
FL0507 
FL0508 
FL0509 
FL0510 
FL0511 
FL0512 
FL0513 
FL0514 
FL0515 
FL0516 
FL0517 
FL0518 
FL0519 
FL0520 
FL0521 
FL0522 
FL0523 
FL0524 
FL0525 
FL0526 
FL0527 
FL0528 
FL0529 
FL0530 
FL0531 
FL0532 
FL0533 
FL0534 
FL0535 
FL0536 
FL0537 
FL0538 
FL0539 
FL0540 
FL0541 
FL0542 
FL0543 
FL0544 
FL0545 
FL0546 
FL0547 
FL0548 
FL0549 
FL0550 
FLQ^53, 


SUBROUTINE  SAFAE 

1  (L,X,U,P,RO,G,E,DU, 

2  US,PS,UIDOT,PIDOT, 
x             FIMZ,ZMDOT, 

3  TENA,FIRO,FIM,FIE, 
IMPLICIT  REAL*8(A-H,0-Z,$) 
DIMENSION  X(L),U(L),P(L),RO(L),G 

1  DG(L),DXSI(L),MIN(L), 

2  US(L),PS(L),UIDOT(L),P 

3  ,TENA(L),FIRO(L),FIM(L) 

4  ,GIP(L),VOL(L),Z(L),DZ< 

5  ,FIMZ(L),ZMDOT(L) 
COMMON  /AB/A(50) 

EQUIVALENCE    ( LL , A( 2) ) , ( T, A( 3) ) , ( 
COMMON   /GAM/GAMA,NG,MU2,G1,G2,G3 

,G12,G13,G14,G15,G16,G 
NG 
MU2 

DETO/QDET,TC,RATE,PCJDET, 
DIFFUS/U2,P2,R02,ARW 
xx*xxxxxxxxxxxxxxxxxxxxxx 


REALX8 
REALX8 
COMMON/ 
COMMON/ 

CXXXXXXXXXXXX 

TN  =  T 
C 

C   INFLON  B.C 
C 


DP,DRO,DG,DXSI,MIN, 

GIP,VOL,Z,DZ) 

(L),E(L),DU(L),DP(L),DRO(L) 

IDOT(L) 
,FIE(L) 
L) 


DT,A(4)),(NCYC,A(12)) 

,G4,G5,G6,G7,G8,G9,G10,G11 

17,G18,G19,G20 


RCJDET,UCJDET,DCJDET,P0DET,RO0DET 
xxxxxxxxxxxxxxxxxxxxx 


.  AT  1=2 


XXXXXXXXXXXXXH 

H 
H 
H 

H 


FL0552 
FL0553 
FL0554 
FL0555 
FL0556 
FL0557 
FL0558 
FL0559 
FL0560 
FL0561 
FL0562 
FL0563 
FL0564 
FL0565 
FL0566 
FL0567 
FL0568 
FL0569 
FL0570 
FL0571 
FL0572 
FL0573 
FL0574 
FL0575 
FL0576 
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U(1)=U2 

P(1)=P2 

R0(1)=R02 

G(1)  =  DSQRT(GAMAXP(1)XR0(D) 

Z(1)=1.D0 

DU(1)=0. 

DP(1)=0. 

DG(1)=0. 

DRO(1)=0. 

c 
c 

DXSI(1)=DXSI(2) 

OUTFLOW    B.C.    AT    I=L 

c 

HFL       SCRIPT    Al 

HFL0577 

HFL0578 

HFL0579 

HFL0580 

HFL0581 

HFL0582 

HFL0583 

HFL0584 

HFL0585 

HFL0586 

HFL0587 

HFL0588 

HFL0589 

U(L)=U(LL)  HFL0590 

P(L)=P(LL)  HFL0591 

RO(L)=RO(LL)  HFL0592 

G(L)=DSQRT(GAMA*P(L)*RO(L))  HFL  0593 

DU(L)=0.  HFL0594 

DPCL)=0.  HFL0595 

DG(L)=0.  HFL0596 

DRO(L)=0.  HFL0597 

DXSI(L)=DXSI(LL)  HFL0598 

C  HFL0599 

RETURN  HFL0600 

EJjjQ HFI  Q£jQJ 

SUBROUTINE  BDOK1 ( L , X, V, DV, MIN, NV)  HFL0602 

IMPLICIT  REAL*S(A-H,0-Z,$)  HFL0603 

DIMENSION  X(L),V(L),DV(L),MIN(L)  HFL0604 

COMMON  /AB/AC50)  HFL0605 

EQUIVALENCE  ( LL , AC  2) ) , ( KEYMON, A( 11 ) )  HFL0606 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX        HFL 06 07 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA  HFL0608 

COMMON  /M0NIT/NC14(4),CASEAV(4),NF16(6),  HFL0609 

1  NM0NU(4),NM0NP(4),NM0NG(4),NM0NR0(4),NM0NZ(4)      HFL0610 

DIMENSION  NM0NV(4,5)  HFL0611 

EQUIVALENCE  (  NMONV(  1 , 1 )  ,  NMONUd  ) )  HFL0612 

C      DIMENSION  NAMEV(5)  HFL0613 

C     DATA  NAMEV/1HU,1HP,1HG,2HR0,1HZ/  HFL0614 

DATA  EPS/1. D-15/  HFL0615 

CXXXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXHFL0616 

C      NV=0  HFL0617 

C      DO  10  N  =  l,5  HFL0618 

C      IF  (NAME.EQ.NAMEV(N))NV=N  HFL0619 

CIO    CONTINUE  HFL0620 

GO  TO  (1,2,3,4,5),  NV  HFL0621 

1  AMIDA=(UMAX-UMIN)**2  HFL0622 
GO  TO  9  HFL0623 

2  AMIDA=(PMAX-PMIN)*X2  HFL0624 
GO  TO  9  HFL0625 

3  AMIDA=((UMAX-UMIN)X(R0MAX-R0MIN))XX2  HFL0626 
GO  TO  9  HFL0627 

4  AMIDA  =  (R0MAX-R0MIN)XX2  HFL0628 
GO  TO  9  HFL0629 

5  AMIDA=1.D0  HFL0630 
GO  TO  9  HFL0631 

9     CONTINUE  HFL0632 

AMIDA=AMIDAXEPSXX2  HFL0633 

AMIDA=EPS  HFL0634 

DO  29  1=2, LL  HFL0635 

ICAT=0  HFL0636 

VLEFT=V(I)-0.5D0XDV(I)  HFL0637 

VRIGHT=V(I)+0.5D0XDV(I)  HFL0638 

VM=V(I-1)  HFL0639 

VP=V(I+1)  HFL0640 

SIGN=(VP-V(I))X(V(I)-VM)  HFL0641 

IF(SIGN.GT.-AMIDA)  GO  TO  22  HFL0642 

21  DV(I)=0.  HFL0643 
ICAT=1  HFL0644 
GO  TO  20  HFL0645 

22  CONTINUE  HFL0646 
SIGN=(VP-VM)XDV(I)  HFL0647 
IF(SIGN.GT.-AMIDA)  GO  TO  24  HFL0648 
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23  DV(I)=0.5D0*(VP-VM) 
VLEFT=V(I)-0.5D0*DV(I) 
VRIGHT=V(I)+0.5D0XDV(I) 
ICAT=2 

24  SIGN=(VLEFT-VM)XDV(I) 
IF(SIGN.GT.-AMIDA)  GO  TO  26 

25  VLEFT=VM 
VRIGHT=2.D0XV(I)-VLEFT 
DV(I)=VRIGHT-VLEFT 
ICAT=3 

26  SIGN=(VP-VRIGHT)*DV(I) 

IF  (SIGN.GT.-AMIDA)  GO  TO  28 

27  VRIGHT=VP 
VLEFT=2.D0XV(I)-VRIGHT 
DV(I)=VRIGHT-VLEFT 
ICAT=3 

28  IFCDABS(DVCD)  . LE . 0 . 5D0XDABSC VP-VM) )    GO    TO    31 

30  DV(I)=0.5D0x(VP-VM) 
ICAT=4 

31  CONTINUE 
20    CONTINUE 

IF  (DV(I)*X2.GT.AMIDA)  GO  TO  40 

DV(I)=0. 
40    CONTINUE 
C      IF  (ICAT.GT.O)  NM0NV(ICAT,NV)=NM0NV(ICAT,NV)+1 
C      IBYTE0=5 
C      MIN(I)=MIN(I) . OR.SHIFT(ICAT,(NV+IBYTEO-2)X3) 

29  CONTINUE 
RETURN 
END 


HFL0649 
HFL0650 
HFL0651 
HFL0652 
HFL0653 
HFL0654 
HFL0655 
HFL0656 
HFL0657 
HFL0658 
HFL0659 
HFL0660 
HFL0661 
HFL0662 
HFL0663 
HFL0664 
HFL0665 
HFL0666 
HFL0667 
HFL0668 
HFL0669 
HFL0670 
HFL0671 
HFL0672 
HFL0673 
HFL0674 
HFL0675 
HFL0676 
HFL0677 
HFLQ678 


SUBROUT 
IMPLICI 

DIMENSI 

COMMON 

EQUIVAL 

Cxxxxxxxx*xxx 
DO  1  1  = 
IM=I-1 
IP=I+1 
DV(I)=0 
1  CONTINU 
RETURN 

EUJtt 


INE  DCOLE(L,X,V,DV,MIN,NV) 

T  REALX8(A-H,0-Z,$) 

ON  X(L),V(L),DV(L),MIN(L) 

/AB/AC50) 

ENCE  (LL,A(2)) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

2,LL 


5X(V(IP)-V(IM)) 


HFL0679 
HFL0680 
HFL0681 
HFL0682 
HFL0683 
XXXXXXXHFL0684 
HFL0685 
HFL0686 
HFL0687 
HFL0688 
HFL0689 
HFL0690 
HFL0691 


SUBROUTINE  PRINT  HFL0692 

(L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  HFL 06 93 

US,PS,UIDOT,PIDOT,  HFL0694 

FIMZ,ZMDOT,  HFL0695 

TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ)  HFL 06 96 

IMPLICIT  REALX8(A-H,0-Z,$)  HFL0697 

DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L),  HFL 06 98 

DG(L),DXSI(L),MIN(L),  HFL0699 

US(L),PS(L),UIDOT(L),PIDOT(L)  HFL 07 0  0 

,TENA(L),FIRO(L),FIM(L),FIE(L)  HFL 07 01 

,GIP(L),VOL(L),Z(L),DZ(L)  HFL 07 02 

,FIMZ(L),ZMDOT(L)  HFL0703 

COMMON  /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT  HFL0704 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  HFL0705 

RSTARL,RSTARR,GSTARL,GSTARR,WL,WR,  HFL0706 

»  CL,CR,CSTARL,CSTARR,UW(6)  HFL0707 

J  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  HFL0708 

LOGICAL  HELEML,HELEMR  HFL0709 

COMMON  /AB/A(50)  HFL0710 

EQUIVALENCE  ( LL , AC  2) ) , (T, A( 3) ) , (NCYC, A( 12) ) , ( DT, A( 4) )  HFL0711 

EQUIVALENCE  (UGAL,A(15))  HFL0712 

COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET  HFL 07 13 

COMMON/ DIFFUS/U2,P2,R02,ARW  HFL 07 14 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 07 15 

,G12,G13,G14,G15,G16,G17,G18,G19,G20  HFL0716 

REALX8  NG  HFL0717 

REALX8  MU2  HFL0718 

COMMON  /M0NIT/NC14(4),CASEAV(4),NF16(6),  HFL0719 

NM0NU(4),NM0NP(4),NM0NG(4),NM0NR0(4),NM0NZ(4)  HFL 07  2  0 


31 


E:  HFL 


SCRIPT 


Al 


DIMENSION  CASAVK4) 

LOGICAL  FULLPR 

COMMON  /PRNT/ROX,PX,UX,ZX 

CX*XXXXX*XXXXXXXXXX*XXXXXX*XX**XXXXXXXX**XXXXXXXX*XXX*XXXXXX*XXXXXXXX 

FULLPR=.TRUE. 
PRINT    1 

1  FORMAT(lHl) 

PRINT    2,    T,DT,NCYC 

2  FORMATdX, 10X, 'RESULTS    AT      T= ', Dll . 5, 5X, ' DT= « , Dll . 5, 5X, ' NCYC= ' , 
1  15//) 

PRINT    3,    AMTOT,ETOT,EKTOT\ 

3  FORMATdX,  • AMTOT= ' , D20 . 14, 


W 


11 


12 


EPTOT,TENTOT 

2X,  'ETOT,EKTOT,EPTOT: 


',3D22.14/ 


IX, 'TENTO 


222 


13 

131 

10 


FORMATdX, 


FORMATdX,  ' 


=  i 


D21 

X 

RO 

DU 

DG 


ZMDOT 

AMDOTN 

AMACH 


14//) 


U 
G 
DP 

DZ1) 
US 

FIMZ 
TEMP 
ENTRO 


DRO 


PS 
AMDOT 
ENTALP 


FORMATdX) 

IF  (UGAL.NE.O.)  PRINT  6,  UGAL 
FORMATC/11X, 'INITIAL  VELOCITY 
DO  10  1=1, L 

10) .NE.l)  GO  TO 


CORRESPONDS  TO   UGAL= ■ , D15 . 6/ ) 


IF  (MOD(I 
PRINT  5 
PRINT  4 
PRINT  44 
PRINT  5 
CONTINUE 
PRINT  12, 


11 


RO(I),G(I),Z(I),DU(I),DP(I),DRO(I) 


X(I),U(I),P(I) 
1  DG(I),DZ(I) 

F0RMATdX,I3,6D12.5,5D11.4) 

ENTRO=P(I)/RO(I)**GAMA 

IF( .NOT. FULLPR)  GO  TO  131 

IF(I.EQ.l)  GO  TO  131 

IM=I-1 

UL  =  UdM)  +  0.5XDUdM) 

PL  =  PdM)  +  0.5XDPdM) 

ROL=ROdM)  +  0.5XDROdM) 

GL=GdM)  +  0.5XDGdM) 

CL=GL/ROL 

ZL=ZdM)  +  0.5XDZdM) 

UR  =  Ud)-0.5xDUd) 

PR=P(I)-0.5*DP(I) 

GR  =  Gd)-0.5XDGd) 

ROR=RO(I)-0.5*DRO(I) 

CR=GR/ROR 

ZR=Z(I)-0.5XDZ(I) 

CALL  RIEMAN(L,I,MIN) 

CALL  FLUXE1(L,I,MIN) 

XI=X(I) 

IFd.EQ.L)  GO 

UX  =  Ud) 

PX  =  Pd) 

ROX=RO(I) 

ZX=Z(I) 

IP=I+1 

XI=0.5D0*(X(I)+X(IP)) 

CONTINUE 

AMACH=UX/DSQRT(GAMAXpX/ROX) 

AMDOT =R0XXUXXCR0SS( XI) 

IF(I.EQ.2)AMDOT0=AMDOT 

AMDOTN=AMDOT/AMD0T0 

ENTALP=(GAMA/(GAMA-1.DO))XPX/ROX+0 

TEMP=PX/(ROXXARW) 

PRINT  13,USd),PSd), 
1  ZMDOT ( I ) ,  FIMZd  ) , AMDOT , AMDOTN , TEMP, ENTALP , AMACH , ENTRO 

F0RMAT(4X,12X,5D12.5,6D11.4) 

CONTINUE 

CONTINUE 


TO  222 


5D0XUXXX2 


HFL0721 
HFL0722 
HFL0723 
XXXHFL0724 
HFL0725 
HFL0726 
HFL0727 
HFL0728 
HFL0729 
HFL0730 
HFL0731 
HFL0732 
HFL0733 
HFL0734 
HFL0735 
HFL0736 
HFL0737 
HFL0738 
HFL0739 
HFL0740 
HFL0741 
HFL0742 
HFL0743 
HFL0744 
HFL0745 
HFL0746 
HFL0747 
HFL0748 
HFL0749 
HFL0750 
HFL0751 
HFL0752 
HFL0753 
HFL0754 
HFL0755 
HFL0756 
HFL0757 
HFL0758 
HFL0759 
HFL0760 
HFL0761 
HFL0762 
HFL0763 
HFL0764 
HFL0765 
HFL0766 
HFL0767 
HFL0768 
HFL0769 
HFL0770 
HFL0771 
HFL0772 
HFL0773 
HFL0774 
HFL0775 
HFL0776 
HFL0777 
HFL0778 
HFL0779 
HFL0780 
HFL0781 
HFL0782 
HFL0783 
HFL0784 
HFL0785 
HFL0786 
HFL0787 
HFL0788 
HFL0789 
HFL0790 
HFL0791 
HFL0792 
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FILE:  HFL 


SCRIPT 


Al 


C   JO 

C 

C 

C 

C40 

C 

C30 

C 

C31 

c 
c 

C301 
C 
C 
C32 

c 
c 
c 

C33 
C 
C 
C 

c 
c 


B  STATI 
DO  40 
CASAV1 
IF  (NC 
CONTIN 
PRINT 
FORMAT 
PRINT 
FORMAT 

1 
PRINT 
FORMAT 

1 
PRINT 
FORMAT 
ICAT0= 
PRINT 

1 
FORMAT 

1 

1 

1 

1 

1 
RETURN 


NE.O)  CASAV1(I)=CASEAV(I)/NC14(I) 


STICS 

1  =  1,4 

U)=0 

1-4(1) 

UE 

30 

(///1X,10(1HX),3X, 'JOB 

31,(NC14(I),I=1,4) 

(IX, 'NO.  OF  VARIOUS  CASES 

4110) 
301,  (CASAVHI),  1  =  1,4) 
(/1X, 'AVERAGE  NUMBER  OF  ITERATIONS  IN 

IX, '   CASAV1(NCASE)=',4(F6.2,4X)) 
32, (NF16(I), 1=1,6) 

(/1X,'N0.  OF  VARIOUS  FLUX  CASES    NF16 ( NFLUX) = • , 61 10 ) 
4 

33,(NMONU(I),I=l,ICAT0),(NMONP(I),I=l,ICAT0), 
(NMONG(I),I=l,ICAT0), ( NMONRO( I ) , 1=1 , ICATO ) 
(/1X,'N0.  OF  MONOTONICITY  INTERVENTIONS  FOR  EACH  VAR.' 

IX, 'IN  EACH  CATEGORY. '/ 

IX, 'NMONU  (ICAT)=»,4I10/ 

IX, 'NMONP  (ICAT)=»,4I10/ 

IX, 'NMONG  (ICAT)=',4I10/ 

IX, 'NMONRO(ICAT)=',4I10/) 


STATISTICS f,3X,10(lHX)//) 

IN  RIEMAN  SOLVER    NC14(NCASE): 

RIEMAN  SOLVER', 


HFL0793 
HFL0794 
HFL0795 
HFL0796 
HFL0797 
HFL0793 
HFL0799 
HFL0800 
HFL0801 
HFL0302 
HFL0803 
HFL0804 
HFL0805 
HFL0806 
HFL0807 
HFL0808 
HFL0809 
HFL0810 
HFL0811 
HFL0812 
HFL0813 
HFL0814 
HFL0S15 
HFL0816 
HFL0817 
HFL0818 


SUBROUTINE  SOF(ISTOP) 

IMPLICIT  REAL*8(A-H,0-Z,$) 

PRINT  1,  ISTOP 

FORMAT ( IX, 3HXXX,3X, ' SOF' , 16 , 3X, 3HXXX/ ) 

XX=-1.D0 

YY=DSQRT(XX) 

CALL  SYSTEM(51,16H  SUBROUTINE  TREE) 

STOP 

END 


HFL0819 
HFL0820 
HFL0821 
HFL0822 
HFL0823 
HFL0824 
HFL0825 
HFL0826 
HFL0827 


Cxxxx 

Cxxxx 

7211 
7201 


SUBROUTINE  RI EMAN( L , I , MIN) 

IMPLICIT  REALX8(A-H,0-Z,$) 

DIMENSION  MIN(L) 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1  RSTARL,RSTARR,GSTARL,GSTARR,WL,WR, 

2  CL,CR,CSTARL,CSTARR,UW(6) 

3  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
LOGICAL  HELEML,HELEMR 

COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR 

2  ,ASTARL,ASTARR 

3  ,RAT,SH 

COMMON  /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 
1  ,XMIN,XMAX,SMIN,SMAX,IVERSA 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11 
1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20 

REALX8  NG 

REALX8  MU2 

COMMON  /AB/A(50) 

COMMON  /M0NIT/NC14(4),CASEAV(4),NF16(6), 
1  NM0NU(4),NM0NP(4),NM0NG(4),NM0NR0(4),NM0NZ(4) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*xxxxxxxxxxxxxxxxxx 

DATA  NMAX/63/ 
DATA  EPS/1. D-8/ 
x*xxxxx*xxxxxxxxxxx*xx**xxxxx*xx*xxxxxxxxx*xxxx*xxxxxxxxxxxxxxxxxx 

IF(PL.GT.0. .AND.PR.GT.O. )  GO  TO  7201 


PRINT  7211, I, PL, PR 

FORMAT( '  FROM  RIEMAN 

CALL  S0F(7654) 

CONTINUE 

UW(6)=1.D20 

WL  =  0. 

WR  =  0. 

ZETAL=PLXXG8 

ZETAR=PRXXG8 

CLG=CL/GAMA 

CRG=CR/GAMA 

ZSTARL=ZL 


I,PL,PR=',I6,2D16.6/) 


HFL0828 
HFL0829 
HFL0830 
HFL0831 
HFL0832 
HFL0833 
HFL0834 
HFL0835 
HFL0836 
HFL0837 
HFL0838 
HFL0839 
HFL0840 
HFL0841 
HFL0842 
HFL0843 
HFL0844 
HFL0845 
HFL0846 
HFL0847 

XHFL0848 
HFL0849 
HFL0850 

XHFL0851 
HFL0852 
HFL0853 
HFL0854 
HFL0855 
HFL0856 
HFL0857 
HFL0858 
HFL0859 
HFL0860 
HFL0861 
HFL0862 
HFL0863 
HFL0864 
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HFL 


SCRIPT 


Al 


TO  102 


+G6XEVERR) 


ZSTARR=ZR 

IF  (ZETAL.LT.ZETAR)  GO 
C   LEFT  PRESSURE  IS  HIGHER 

101  CONTINUE 
EVERR=(PL-PR)/PR 
USR=UR+CRG*EVERR/DSQRT(1 
SR=USR 

UEL=UL-G7XCL*CZETAR-ZETAL)/ZETAL 
SL=UEL 

NL=2 
NR=2 

IF  (USR.GE.UL)  NL=1 
IF  (UEL.LE.UR)  NR=1 
IF  (DABS(EVERR) .LT.EPS)  GO  TO  100 
IF  (NL.EQ.2.AND.NR.EQ.1)  GO  TO  7001 
GO  TO  100 
C   RIGHT  PRESSURE  IS  HIGHER 

102  CONTINUE 
EVERL=(PR-PL)/PL 
USL=UL-CLG*EVERL/DSQRT(1 
SL=USL 

UER=UR+G7*CRX(ZETAL-ZETAR)/ZETAR 

SR=UER 

NL=2 

NR  =  2 

IF  (UER.GE.UL)  NL=1 

(USL.LE.UR)  NR=1 

(DABS(EVERL) .LT.EPS)  GO  TO  100 


+G6XEVERL) 


100 


IF 
IF 
IF 
GO 


(NL.EQ 
TO  100 
CONTINUE 
IF  (NL.EQ 
IF  (NL.EQ 
IF  (NL.EQ 
IF  (NL.EQ 


1.AND.NR.EQ.2)  GO  TO  7001 


NCASE=1 
NCASE=2 
NCASE=3 
NCASE=4 

LT.EPSX(PMAX-UMIN)) 


1.AND.NR.EQ.2) 

2.AND.NR.EQ.2) 

2.AND.NR.EQ.1) 

l.AND.NR.EQ.l) 
I  F(  DABS  (  PL-PR)  +  DABS(UL-UR).LT.EPSX(PMAX-UMIN))    NCASE=4 
UMIDA  =  EPSxDMAXHCL,CR) 
DUDZL=-G7*CL/ZETAL 
DUDZR=    G7XCR/ZETAR 

ZETA=(-(UR-UL)+ZETARXDUDZR-ZETALXDUDZL)/(DUDZR-DUDZL) 
IF    (ZETA.LE.O. )    GO   TO    7002 
N  =  0 

GO  TO  (1,2,3,4),  NCASE 
:   THE  CASE   ES 

I  ITYPE=NCASE 
HELEML=. FALSE. 
HELEMR=.TRUE. 

II  N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 

ZETAF=ZETA 

UEL=UL-G7XCLX(ZETAF-ZETAL)/ZETAL 

PPR=(ZETAF/ZETAR)XXNG 

EVERR=PPR-1. 

SQRR=DSQRT(1.+G6*EVERR) 

USR=UR+CRGXEVERR/SQRR 

DU=UEL-USR 

IF  (DABS(DU) .LE.UMIDA)  GO  TO  10 

DUDZR=NG*CRGX(PPR/ZETAF)*(1.+G9*EVERR)/SQRRXX3 
ZETA=ZETAF+DU/(DUDZR-DUDZL) 
GO  TO  11 
10    CONTINUE 

USTAR=(UEL+USR)/2. 
PSTAR=PPRXpR 
CSTARL=CL+(UL-USTAR)/G7 
RSTARL=GAMAXPSTAR/CSTARLXX2 
GSTARL=CSTARLXRSTARL 
:   EQU.  NO.   69.01  OF  THE  BOOK  BY  COURANT-FRIEDRICHS . 
WWR=G11X(USTAR-UR)XR0R 
WR=HWR+DSQRT(GRXX2+WWRXX2) 
RSTARR=RORxWR/(WR-RORX(USTAR-UR)) 
GSTARR=DSQRT(GAMAXPSTARXRSTARR) 
CSTARR=GSTARR/RSTARR 


HFL0865 
HFL0866 
HFL0867 
HFL0868 
HFL0869 
HFL0870 
HFL0871 
HFL0872 
HFL0873 
HFL0874 
HFL0875 
HFL0876 
HFL0877 
HFL0878 
HFL0879 
HFL0880 
HFL0881 
HFL0882 
HFL0883 
HFL0884 
HFL0885 
HFL0886 
HFL0887 
HFL0888 
HFL0889 
HFL0890 
HFL0891 
HFL0892 
HFL0893 
HFL0894 
HFL0895 
HFL0896 
HFL0897 
HFL0898 
HFL0899 
HFL0900 
HFL0901 
HFL0902 
HFL0903 
HFL0904 
HFL0905 
HFL0906 
HFL0907 
HFL0908 
HFL0909 
HFL0910 
HFL0911 
HFL0912 
HFL0913 
HFL0914 
HFL0915 
HFL0916 
HFL0917 
HFL0918 
HFL0919 
HFL0920 
HFL0921 
HFL0922 
HFL0923 
HFL0924 
HFL0925 
HFL0926 
HFL0927 
HFL0928 
HFL0929 
HFL0930 
HFL0931 
HFL0932 
HFL0933 
HFL0934 
HFL0935 
HFL0936 
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ILE:  HFL       SCRIPT    Al 

WRE=WR/ROR+UR  HFL0937 

UN(1)=UL-CL  HFL0938 

UW(2)=USTAR-CSTARL  HFL0939 

UW(3)=USTAR  HFL0940 

UW(4)=WRE  HFL0941 

UW(5)=HRE  HFL0942 

GO  TO  5  HFL0943 

C   THE  CASE   SS  HFL0944 

2  ITYPE=NCASE  HFL0945 
HELEML=.TRUE.  HFL0946 
HELEMR=.TRUE.  HFL0947 

21    N=N+1  HFL0948 

IF  (N.GT.NMAX)  GO  TO  7003  HFL0949 

ZETAF=ZETA  HFL0950 

PF=ZETAF**NG  HFL0951 

PPL-PF/PL  HFL0952 

PPR=PF/PR  HFL0953 

EVERL=PPL-1.  HFL0954 

EVERR  =  PPR-1  .  HFL0955 

SQRL=DSQRT(1.+G6*EVERL)  HFL0956 

SQRR=DSQRT(1.+G6*EVERR)  HFL0957 

USL=UL-CLG*EVERL/SQRL  HFL0958 

USR=UR+CRG*EVERR/SQRR  HFL0959 

DU=USL-USR  HFL0960 

IF  (DABS(DU) .LE.UMIDA)  GO  TO  20  HFL0961 

DUDZL=-NG*CLG*(PPL/ZETAF)*(1 . +G9*EVERL)/SQRL*X3  HFL 0962 

DUDZR=  NG*CRGX(PPR/ZETAF)X(1 . +G9*EVERR)/SQRRXX3  HFL0963 

ZETA=ZETAF+DU/(DUDZR-DUDZL)  HFL 096 4 

GO  TO  21  HFL0965 

20    CONTINUE  HFL0966 

USTAR=(USL+USR)/2.  HFL0967 

PSTAR=(PPL*PL+PPR*PR)/2.  HFL0968 

WWR=G11*(USTAR-UR)XR0R  HFL0969 

WR=WWR+DSQRT(GRXX2+WWRXX2)  HFL 097  0 

WWL=-Gllx(USTAR-UL)xROL  HFL0971 

WL=WWL+DSQRT(GL*X2+WWL*X2)  HFL 0  97  2 

RSTARL=ROLXWL/(WL  +  ROLX(USTAR-UD)  HFL  097  3 

RSTARR=ROR*WR/(WR-RORX(USTAR-UR))  HFL 0  97 4 

GSTARL=DSQRT(GAMAXPSTARXRSTARL)  HFL 0  97 5 

GSTARR=DSQRT(GAMAXPSTARXRSTARR)  HFL 0  976 

CSTARL=GSTARL/RSTARL  HFL0977 

CSTARR=GSTARR/RSTARR  HFL0978 

WLE=-HL/ROL+UL  HFL0979 

WRE=WR/ROR+UR  HFL0980 

UW(1)=WLE  HFL0981 

UW(2)=HLE  HFL0982 

UW(3)=USTAR  HFL0983 

UW(<4)=WRE  HFL0984 

UW(5)=WRE  HFL0985 

GO  TO  5  HFL0986 

C   THE  CASE   SE  HFL0987 

3  ITYPE=NCASE  HFL0988 
HELEML=.TRUE.  HFL0989 
HELEMR=. FALSE.  HFL0990 

31    N=N+1  HFL0991 

IF  (N.GT.NMAX)  GO  TO  7003  HFL0992 

ZETAF=ZETA  HFL0993 

UER=UR+G7XCRX(ZETAF-ZETAR)/ZETAR  HFL 0994 

PPL=(ZETAF/ZETAL)*XNG  HFL0995 

EVERL=PPL-1.  HFL0996 

SQRL=DSQRT(1.+G6XEVERL)  HFL0997 

USL=UL-CLG*EVERL/SQRL  HFL0998 

DU=USL-UER  HFL0999 

IF  (DABSCDU) .LE.UMIDA)  GO  TO  30  HFL1000 

DUDZL=-NG*CLGX(PPL/ZETAF)X(1.+G9XEVERL)/SQRL*X3  HFL 10 01 

ZETA=ZETAF+DU/(DUDZR-DUDZL)  HFL 1002 

GO  TO  31  HFL1003 

30    CONTINUE  HFL1004 

USTAR=(USL+UER)/2.  HFL1005 

PSTAR=PPLXPL  HFL1006 

CSTARR=CR-(UR-USTAR)/G7  HFL1007 

RSTARR=GAMA*PSTAR/CSTARRXX2  HFL  10  08 
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HFL 


SCRIPT 


Al 


GSTARR=CSTARR*RSTARR  HFL1009 

WWL=-G11*(USTAR-UL)XR0L  HFL1010 

WL=WWL  +  DSQRT.(GLX*2+WWLx*2)  HFL  1011 

WLE=-WL/ROL+Ul  HFL1012 

RSTARL=ROLXWL/(WL+ROLX(USTAR-UL))  HFL  101 3 

GSTARL=DSQRT(GAMAXPSTARXRSTARL)  HFL 1014 

CSTARL=GSTARL/RSTARL  HFL1015 

UW(1)=WLE  HFL1016 

UW(2)=WLE  HFL1017 

UW(3)=USTAR  HFL1018 

UW(4)=USTAR+CSTARR  HFL1019 

UW(5)=UR+CR  HFL1020 

GO  TO  5  HFL1021 

C   THE  CASE   EE  HFL1022 

4  ITYPE=NCASE  HFL1023 
HELEML=. FALSE.  HFL1024 
HELEMR=. FALSE.  HFL1025 
PSTAR=ZETA*XNG  HFL1026 
USTAR=UL-G7XCLX(ZETA-ZETAL)/ZETAL  HFL 1027 
CSTARL=CL+(UL-USTAR)/G7  HFL1028 
CSTARR=CR-(UR-USTAR)/G7  HFL1029 
RSTARL=GAMAXPSTAR/CSTARLXX2  HFL 10  30 
RSTARR=GAMAXPSTAR/CSTARRXX2  HFL 1031 
GSTARL=RSTARLXCSTARL  HFL1032 
GSTARR=RSTARRXCSTARR  HFL1033 
UW(1)=UL-CL  HFL1034 
UW(2)=USTAR-CSTARL  HFL1035 
UW(3)=USTAR  HFL1036 
UW(4)=USTAR+CSTARR  HFL1037 
UW(5)=UR+CR  HFL1038 
N=l  HFL1039 
GO  TO  5  HFL1040 

5  CONTINUE  HFL1041 
DO  6  K=l,6  HFL1042 
NFLUX=K  HFL1043 
IF  (UW(K).GE.O.)  GO  TO  61  HFL1044 

6  CONTINUE  HFL1045 
NFLUX=6  HFL1046 

61    CONTINUE  HFL1047 

MIN(I)=MIN(I) .OR.NCASE  HFL1048 

MIN(I)=MIN(I) . OR. SHI  FT (N, 3)  HFL 1049 

MIN(I)=MIN(I) . OR. SHI  FT (NFL UX, 9)  HFL 1050 

NC1<4(NCASE)=NC14(NCASE)  +  1  HFL  10  51 

CASEAV(NCASE)=CASEAV(NCASE)+N  HFL 1052 

NF16(NFLUX)=NF16(NFLUX)+1  HFL 1053 

IF(A(3) .NE.O. )G0  TO  666  HFL1054 

IF(I.NE.2)  GO  TO  666  HFL1055 
PRINT  66  7,I,NFLUX,NCASE,PL,UL,R0L,PR,UR,R0R,USTAR,PSTAR,RSTARL,    HFL 1056 

1  RSTARR,(KK,UW(KK),KK=1,6)  HFL1057 
FORMAT (/ IX, 'I,NFLUX,NCASE=',3I5/1X, • PL , UL , ROL , PR, UR, ROR= ■ , 6D12 . 4/  HFL 1058 

1  IX, 'USTAR,PSTAR,RSTARL,RSTARR=',4D13.4/  HFL1059 

2  IX, 'KK,UW(KK)=»,6(I4,2X,D13.4)/)  HFL1060 
CONTINUE  HFL1061 
RETURN  HFL1062 
CONTINUE  HFL1063 
PRINT  7101,  PL,UL,PR,UR,ZETAL,ZETAR,SL,SR,NL,NR,I  HFL1064 
F0RMAT(//1X, 'FROM  RIEMAN.   AN  IMPOSSIBLE  CASE  OF  EXPANSION/SHOCK'  HFL1065 

1  //1X, »PL,UL,PR,UR=',4D25.14//  HFL1066 

2  IX,  'ZETAL,ZETAR,SL,SR=»,4D25.1<+//  HFL1067 

3  IX, !NL,NR,I=',3I10//)  HFL1068 
CALL  SOFC7001)  HFL1069 
CONTINUE  HFL1070 
PRINT  7102,  ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR,N,NCASE,I  HFL1071 
F0RMAT(//1X, 'FROM  RIEMAN.   NEGATIVE  PRESSURE  AT  THE  INTERSECTION' , HFL1072 

•OF  L  AND  R  EXPANSION  BRANCHES'//  HFL1073 

•IT  MEANS  THAT  A  CAVITATION  TENDS  TO  FORM.   THIS',     HFL1074 

IX, 'POSSIBILITY  IS  EXCLUDED  IN  PRESENT  VERSION'//  HFL1075 

IX, 'ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR=',9DI0.3//   HFL 1076 

IX, 'N,NCASE,I=',3I10//)  HFL1077 

CALL  SOF(7002)  HFL1078 

7003  CONTINUE  HFL1079 

PRINT  7103,  I,N,NCASE,DU,UMIDA,EPS,PL,UL,PR,UR,  HFL1080 


667 

666 
7001 

7101 

7002 
7102 


IX, 
IX, 
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1            ZETA,ZETAF,ZETAL,ZETAR,DUDZL,DUDZR  HFL1081 

7103  FORMATC//1X, 'FROM  RIEMAN.   NUMBER  OF  ITERATIONS  EXCEEDED.1//  HFL1082 

1  IX, ,I,N,NCASE,DU,UMIDA,EPS=,,3I6,3D18.6//  HFL1083 

2  "  IX, 'PL,UL,PR,UR,ZETA,ZETAF=» ,6018.10//  HFL1084 

3  IX, 'ZETAL,ZETAR,DUDZL,DUDZR=',4D18.10//)  HFL1085 
CALL  SOF(7003)  HFL1086 
RETURN  HFL1087 

Eii£ HFI  L0&2 

SUBROUTINE  MAGA( L , I , MIN)  HFL1089 

IMPLICIT  REAL*8(A-H,0-Z,$)  HFL1090 

DIMENSION  MIN(L)  HFL1091 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 1092 

1           ,G12,G13,G14,G15,G16,G17,G18,G19,G20  HFL1093 

REAL*8  NG  HFL1094 

REAL*8  MU2  HFL1095 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  HFL1096 

1  RSTARL,RSTARR,GSTARL,GSTARR,WL,WR,  HFL1097 

2  CL,CR,CSTARL,CSTARR,UW(6)  HFL1098 

3  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  HFL1099 
LOGICAL  HELEML,HELEMR  HFL1100 
COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR  HFL1101 

2  ,ASTARL,ASTARR  HFL1102 

3  ,RAT,SH  HFL1103 
COMMON  /GRADS/DUDXIL,DPDXIL,DGDXIL,DRDXIL,  HFL1104 

1              DUDXIR,DPDXIR,DGDXIR,DRDXIR  HFL1105 

COMMON  /AB/AC50)  HFL1106 

REAL  LU,LP,LRO  HFL1107 

DATA  EPS/1. D-6/  HFL1108 

C   WE  HERE  SOLVE  FOR  THE  TIME-DERIVATIVES  ALONG  THE  CONTACT  SURFACE,     HFL1110 

C   NAMELY   DUIDT,DPIDT.   FROM  THESE  WE  ALSO  OBTAIN  THE  OTHER  HFL1111 

C   TIME-DERIVATIVES  (SEE  COMMON  /STEP1/).  HFL1112 
C   WE  COMPUTE  THE  COEFFICIENTS  FOR  TWO  EQUATIONS  FOR  DUIDT,DPIDT.   THESEHFL1113 

C   ARE         AAL*DUIDT+BBL*DPIDT=CCL  HFL1114 

C              AARXDUIDT+BBR*DPIDT=CCR  HFL1115 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXHF 11116 

C     IF(SH.LE.EPS)RAT=0.  HFL1117 

IF  C.NOT.HELEML)  GO  TO  12  HFL1118 

11  CONTINUE  HFL1119 
DP=PSTAR-PL  HFL1120 
IF  (DABS(DP) .LT.EPSXPSTAR)  GO  TO  12  HFL1121 
DU=USTAR-UL  HFL1122 
Z1=1./DP  HFL1123 
Z2=0.5/(PSTAR+MU2*PL)  HFL1124 
ZA=Z1-Z2  HFL1125 
ZB=Z1+MU2*Z2  HFL1126 
LU=-ROL*DU*(GAMAXPL*ZB+0.5)+WL  HFL 11 27 
LRO=0.5*DP/ROL  HFL1128 
LP=1.+ZB*DP  HFL1129 
AAL=1.+ZA*DP  HFL1130 
3BL=-ZA*DU+WL/GSTARLXX2  HFL1131 
CCL=-LU*DUDXIL-LRO*DRDXIL-LPXDPDXIL  HFL 1132 
CCL=CCL-WL*USTARXRAT/RSTARL  HFL 11 33 

1   +ULXRATXDUXCGAMAXPLXZB+0.5)  HFL1134 

GO  TO  10  HFL1135 

12  CONTINUE  HFL1136 
A1=DUDXIL+DPDXIL/GL  HFL1137 
BETA=GSTARL/GL  HFL1138 
ASTARL=A1+(G3/GL)X(CLXDGDXIL-G4XDPDXIL)X(BETAXXG5-1.)  HFL 11 39 
AAL=1.  HFL1140 
BBL=1./GSTARL  HFL1141 
CCL=-GSTARLXASTARL/DSQRT(BETA)  HFL 11 42 
GE0M=RAT*((GAMA-1. )*UL+2.XCL)X  HFL 11 43 

1  (BETAXXG13-1 . )/(R0LX(GAMA-3. ))  HFL1144 

1   -4.XRATXCLXCBETAXXG14-1 . )/(R0LX(3.XGAMA-5. ))  HFL1145 

ASTARL=ASTARL-GEOM  HFL1146 

EVER1=  GSTARLXGEOM/DSQRTCBETA)  HFL1147 

EVER2=-RAT*USTARXCSTARL  HFL1148 

CCL=CCL+EVER1+EVER2  HFL1149 

GO  TO  10  HFL1150 

10    CONTINUE  "E!"Hfi 

IF  ( .NOT.HELEMR)  GO  TO  22  HFL1152 
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21  CONTINUE  HFL1153 
DP  =  PSTAR-PR  HFL1154 
IF  CDABS(DP).LT.EPSXPSTAR)  GO  TO  22  HFL1155 
DU=USTAR-UR  HFL1156 
Z1=1./DP  HFL1157 
Z2=0.5/(PSTAR+MU2*PR)  HFL1158 
ZA=Z1-Z2  HFL1159 
ZB=Z1+MU2XZ2  HFL1160 

LU=-ROR*DU*(GAMA*PRXZB+0.5)-WR  HFL 1161 

LRO=0.5*DP/ROR  HFL1162 

LP=1.+ZB*DP  HFL1163 

AAR=1.+ZA*DP  HFL1164 

BBR=-ZA*DU-WR/GSTARRX*2  HFL1165 

CCR=-LU*DUDXIR-LROXDRDXIR-LPxDPDXIR  HFL 1166 

CCR=CCR+WR*USTAR*RAT/RSTARR  HFL 1167 

1   +URxRAT*DUX(GAMAxpR*ZB+0.5)  HFL1168 

GO  TO  20  HFL1169 

22  CONTINUE  HFL1170 
A1=DUDXIR-DPDXIR/GR  HFL1171 
BETA=GSTARR/GR  HFL1172 
ASTARR=A1-(G3/GR)X(CRXDGDXIR-G4*DPDXIR)X(BETA**G5-1 .  )  HFL1173 
AAR=1.  HFL1174 
BBR=-1./GSTARR  HFL1175 
CCR=GSTARR*ASTARR/DSQRT(BETA)  HFL 1176 
GE0M=RATX(-(GAMA-1.)XUR+2.*CR)X(BETA**G13-1.)  HFL 1177 

1  /(R0RX(GAMA-3.))  HFL1178 

2  -4.XRATxCRX(BETAXXG14-l.)/CR0RX(3.XGAMA-5.))  HFL1179 
ASTARR=ASTARR+GEOM  HFL1180 
EVER1=GSTARRXGE0M/DSQRT(BETA)  HFL 1181 

EVER2=RATXUSTAR*CSTARR  HFL1182 

CCR=CCR+EVER1+EVER2  HFL1183 

GO  TO  20  HFL1184 

20    CONTINUE  HFL1185 

DET=AALXBBR-AARXBBL  HFL1186 

DUIDT=(CCLXBBR-CCRXBBL)/DET  HFL 1187 

DPIDT=-(CCLXAAR-CCR*AAL)/DET  HFL 1188 

DRIDTL=DPIDT/CSTARL**2  HFL1189 

DRIDTR=DPIDT/CSTARR**2  HFL1190 

DGIDTL=0.5*GSTARL*(DPIDT/PSTAR+DRIDTL/RSTARL)  HFL 1191 

DGIDTR=0.5XGSTARRX(DPIDT/PSTAR+DRIDTR/RSTARR)  HFL 1192 

RETURN  HFL1193 

END HF1  1194 

SUBROUTINE  FLUXEC L , I , MIN)  HFL1195 

IMPLICIT  REAL*8(A-H,0-Z,$)  HFL1196 

DIMENSION  MIN(L)  HFL1197 

COMMON  /AB/AC50)  HFL1198 

EQUIVALENCE  ( DT, A(4) ) , (NCYC, A( 12) )  HFL1199 

COMMON  /GAM/GAMA,NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11  HFL 120  0 

1  ,G12,G13,G1<+,G15,G16,G17,G18,G19,G20  HFL1201 

REALX8  NG  HFL1202 

REALX8  MU2  HFL1203 

COMMON  /GRADS/DUDXIL,DPDXIL,DGDXIL,DRDXIL,  HFL1204 

1  DUDXIR,DPDXIR,DGDXIR,DRDXIR  HFL1205 

COMMON  /STEP0/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR,  HFL1206 

1  RSTARL,RSTARR,GSTARL,GSTARR,WL,WR,  HFL1207 

2  CL,CR,CSTARL,CSTARR,UW(6)  HFL1208 

3  ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR  HFL1209 
LOGICAL  HELEML,HELEMR  HFL1210 
COMMON  /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR  HFL1211 

2  ,ASTARL,ASTARR  HFL1212 

3  ,RAT,SH  HFL1213 
COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,P0DET,RO0DET  HFL 121 4 
COMMON  /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN  HFL1215 

1  ,GIH  HFL1216 

2  ,FIH4,ZMD0TL,ZMD0TR  HFL1217 
COMMON  /PRNT/ROX,PX,UX,ZX  HFL1218 

ENTRY  FLUXE1(L,I,MIN)  HFL1220 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXHFL1221 

DT2=DT/2.  HFL1222 

GO  TO  (1,2,3,4,5,6),NFLUX  HFL1223 

1     CONTINUE  HFL1224 
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DUDXIX=DUDXIL  HFL1225 

DRDXIX=DRDXIL  HFL1226 

DPDXIX-=PPDXIL  HFL1227 

DGDXIX=0.5*GLX(DPDXIL/PL+DRDXIL/R0n  HFL 1228 

DUDTX=-DPDXIL  HFL1229 

DR0DTX=-R0Lx*2XDUDXIL  HFL1230 

DPDTX=-GLX*2xDUDXIL  HFL1231 

DRODTX=DRODTX-RATXROLXUL  HFL1232 

DPDTX=DR0DTX*CL**2  HFL1233 

DGDTX=G6XGLXDPDTX/PL  HFL1234 

UX  =  UL  HFL1235 

PX=PL  HFL1236 

ROX=ROL  HFL1237 

ZX  =  ZL  HFL1233 

GX  =  GL  HFL1239 

GO  TO  9  HFL1240 

CONTINUE  HFi 1241 

DUDXIX=DUDXIR  HFL1242 

DPDXIX=DPDXIR  HFL1243 

DRDXIX=DRDXIR  HFL1244 

DGDXIX=0.5xGRX(DPDXIR/PR+DRDXIR/ROR)  HFL 1245 

DUDTX=-DPDXIR  HFL1246 

DPDTX=-GRX*2xDUDXIR  HFL1247 

DR0DTX=-R0RX*2XDUDXIR  HFL1248 

DRODTX=DRODTX-RAT*RORXUR  HFL1249 

DPDTX=DR0DTXXCRX*2  HFL1250 

DGDTX=G6XGRXDPDTX/PR  HFL1251 

UX=UR  HFL1252 

PX=PR  HFL1253 

ROX=ROR  HFL1254 

ZX  =  ZR  HFL1255 

GX=GR  HFL1256 

GO  TO  9  HFL1257 

CONTINUE  HFL1253 

BETA0=(MU2X(UL/CL+G7))**(1 ./MU2)  HFL1259 

A1=DUDXIL+DPDXIL/GL  HFL1260 

AO=Al+(G3/GL)x(CLxDGDXIL-G4XDPDXIL)X(BETAOX*G5-l.)  HFL  1261 

EVERl=-(CGAMA-l)*UL+2xCL)x(BETA0*XG13-l. )/(GAMA-3. )  HFL 126 2 

EVER2=4.xCLX(BETA0x*G14-l  .  )/(3 . XGAMA-5 .  )  HFL  126  3 

EVER=(EVER1+EVER2)XRAT/R0L  HFL 126 4 

AO=(AO+EVER)  HFL1265 

DUDAX=AO  HFL1266 

DPDAX=GLXBETAOXAO  HFL1267 

C0=MU2X(UL+G7XCL)  HFL1268 

UX=C0  HFL1269 

ROX=GLXBETAO/CO  HFL1270 

ZX=ZSTARL  HFL1271 

PX=R0XXC0XX2/GAMA  HFL1272 

GX=ROXXCO  HFL1273 

DRODAX=ROLXBETAOXX(1./G4)  HFL 127  4 

1  X((DRDXIL/R0L-DPDXIL/(GAMAXPL))XDSQRT(BETA0)  HFL1275 

2  +AO/CO)  HFL1276 
DGDAX=BETA0XDSQRT(BETA0)X(DGDXIL-G4XDPDXIL/CL)  HFL 1277 

1      +G4XROLXAOXBETAOXX(1./G4)  HFL1278 

DPDAX=DPDAX+RATXUXXCOXDSQRT(BETAO)  HFL 127  9 

G41  =  l ./G4+0.5  HFL1280 

DRODAX=(DRDXIL-DPDXIL/(CLXCD)    XBETA0x*G41+DPDAX/(C0xC0)  HFL  1281 

DGDAX=0.5XGAMAX(PXXDRODAX+ROXXDPDAX)/GX  HFL 128  2 

DUDBX=-CLXBETA0X*(-1./G4)/G4  HFL 128 3 

DPDBX=PLXBETA0XXMU2/G6  HFL1284 

DRODBX=ROLXBETAOXX(-MU2)/G4  HFL 128  5 

DGDBX=GL  HFL1286 

GO  TO  9  HFL1287 

CONTINUE  HFL1288 

BETA0=(MU2X(-UR/CR+G7))XX(1./MU2)  HFL 128  9 

A1=DUDXIR-DPDXIR/GR  HFL1290 

A0=A1+(G3/GR)X(-CRXDGDXIR+G4XDPDXIR)X(BETA0XXG5-1. )  HFL 1291 

EVER1=(-(GAMA-1 . )XUR+2XCR)X( BETA0X*G13-1 . )/(GAMA-3. )  HFL 1292 

EVER2=-4.XCRX(BETA0X*G14-1 . )/( 3 . XGAMA-5 . )  HFL 129 3 

EVER=(EVER1+EVER2)XRAT/R0R  HFL 129  4 

A0=(A0+EVER)  HFL1295 

DUDAX=AO  HFL1296 
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DPDAX=-GRXBETAOXAO  HFL1297 

C0=MU2*(-UR+G7*CR)  HFL1298 

DRODAX=-R0RXBETA0XX(l./G4)  HFL  1299 

1  *a-DRDXIR/ROR+DPDXIR/(GAMA*PR))xDSQRT(BETA0)  HFL1300 

2  +AO/CO)  HFL1301 
UX=-CO  HFL1302 
ROX=GRXBETA0/C0  HFL1303 
ZX=ZSTARR  HFL1304 
PX=ROXxC0XX2/GAMA  HFL1305 
GX=ROXxC0  HFL1306 

DGDAX=BETA0*DSQRT(BETA0)*(-DGDXIR+G4*DPDXIR/CR)  HFL 1307 

1      +G4XRORXA0XBETA0xx(l ./G4)  HFL1308 

DGDAX=-DGDAX  HFL1309 

DPDAX=DPDAX-RATXUXxCOXDSQRT(BETAO)  HFL 1310 

G41=l./G4+0.5  HFL1311 

DRODAX=(DRDXIR-DPDXIR/(CRXCR))    xBETA0xxG41+DPDAX/(C0xC0)  HFL  131 2 

DGDAX=0.5*GAMAX(PX*DRODAX+ROX*DPDAX)/GX  HFL 131 3 

DUDBX=CRXBETA0XX(-1  ./G4VG4  HFL1314 

DPDBX=PRXBETA0XXMU2/G6  HFL1315 

DRODBX=RORXBETA0xx(-MU2)/G4  HFL  1316 

DGDBX=GR  HFL1317 

GO  TO  9  HFL1318 

3  CONTINUE  HFL1319 
UX=USTAR  HFL1320 
PX=PSTAR  HFL1321 
ROX=RSTARL  HFL1322 
ZX=ZSTARL  HFL1323 
GX=GSTARL  HFL1324 
DUDXIX=-DPIDT/GSTARLxx2  HFL1325 
DPDXIX=-DUIDT  HFL1326 
DUDXIX=DUDXIX-RATXUSTAR/RSTARL  HFL 1327 
IF  (.NOT.HELEML)  GO  TO  32  HFL1328 

31  CONTINUE  HFL1329 
DRDXIX=(RSTARL/WL)X*2X(3.XDUIDT+DPIDTX(1.+3.X(WL/GSTARL)**2)/WL    HFL 1330 

1  +DUDXILxWLx((GL/WL)xx2+3.)+3.xDPDXIL  HFL1331 

2  +DRDXILx(WL/R0L)xx2)  HFL1332 
EVER1=ULXRSTARLXX2XRATX((GL/WL)XX2+1.)/(R0L^WL)  HFL 1333 
EVER2=2.XRSTARLXUSTARXRAT/WL  HFL 1334 
DRDXIX=DRDXIX+EVER1+EVER2  HFL 1335 
GO  TO  33  HFL1336 

32  CONTINUE  HFL1337 
BETA=GSTARL/GL  HFL1338 
SQB=DSQRT(BETA)  HFL1339 
DRODA  =  ROLXBETAXX(1./G4)X(SQBX(DRDXIL/ROL-DPDXIL/(GAMAXPD)  HFL  1340 

1                         +ASTARL/CSTARL)  HFL1341 

DRDXIX=DR0DA/SQB+DPIDT/(GSTARLXCSTARLXX2)  HFL 1342 

DPDA=GSTARLX(ASTARL+RATXUSTARXCSTARL/(GLX    5QB))  HFL1343 

G41=l./G4+0.5  HFL1344 

DR0DA  =  (DRDXIL-DPDXIL/(CLXCD)    XBETAXXG41+DPDA/(CSTARLXX2)  HFL  1345 

DRDXIX=    DR0DA/SQB+DPIDT/CGSTARLXCSTARLXX2)  HFL1346 

33  CONTINUE  HFL1347 
DGDXIX=0 .5XGXX(DPDXIX/PX+DRDXIX/R0X)  HFL1348 
DUDTX=DUIDT  HFL1349 
DPDTX=DPIDT  HFL1350 
DGDTX=G6XGSTARLXDPIDT/PSTAR  HFL 1351 
DR0DTX=DPIDT/CSTARLXX2  HFL1352 
GO  TO  9  HFL1353 

4  CONTINUE  HFL1354 
UX=USTAR  HFL1355 
PX=PSTAR  HFL1356 
ROX=RSTARR  HFL1357 
ZX=ZSTARR  HFL1358 
GX=GSTARR  HFL1359 
DUDXIX=-DPIDT/GSTARRXX2  HFL1360 
DUDXIX=DUDXIX-RATXUSTAR/RSTARR  HFL 1361 
DPDXIX=-DUIDT  HFL1362 
IF  ( .NOT.HELEMR)  GO  TO  42  HFL1363 

41  CONTINUE  HFL1364 
DRDXIX=(RSTARR/WR)XX2X(3.XDUIDT-DPIDTX(1 . +3 . x( WR/GSTARR)XX2)/WR    HFL 136 5 

1  -DUDXIRXWRX((GR/WR)xx2+3. )+3.XDPDXIR  HFL1366 

2  +DRDXIRX(WR/R0R)XX2)  HFL1367 
EVER1=URXRSTARRXX2XRATX((GR/WR)XX2+1 . V(RORXWR)  HFL1368 
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ILE:  HFL       SCRIPT    Al 

EVER2=2.*RSTARR*USTAR*RAT/WR  HFL 136 9 

DRDXIX=DRDXIX-EVER1-EVER2  HFL 137  0 

GO  TO  43  HFL1371 

42  CONTINUE  HFL1372 
BETA=GSTARR/GR  HFL1373 
SQB=D5QRT(BETA)  HFL1374 
DR0DA=-R0R*BETAX*(1 ./G4)X( SQB*( -DRDXIR/ROR+DPDXIR/(GAMA*PR) )  HFL 137 5 

1  +ASTARR/CSTARR)  HFL1376 

DRDXIX=DR0DA/SQB-DPIDT/(GSTARR*CSTARRx*2)  HFL 137 7 

DPDA=-GSTARR*(ASTARR+RAT*USTAR*CSTARR/(GRX    SQB))  HFL1373 

G41  =  l ./G4+0 .5  HFL1379 

DRODA=(DRDXIR-DPDXIR/(CR*CR))    *BETA**G41+DPDA/ (CSTARRXX2 )  HFL 133  0 

DRDXIX=    DR0DA/SQB-DPIDT/(GSTARRXCSTARRX*2)  HFL1381 

43  CONTINUE  HFL1382 

DGDXIX=0.5*GX*(DPDXIX/PX+DRDXIX/R0X)  HFL 138  3 

DUDTX=DUIDT  HFL1384 

DPDTX=DPIDT  HFL1385 

DGDTX=G6XGSTARRXDPIDT/PSTAR  HFL 1386 

DR0DTX=DPIDT/CSTARRx*2  HFL1387 

GO  TO  9  HFL1338 

9  CONTINUE  HFL13S9 
RETURN  HFL1390 

CXXXXXXXXXX*XXXXXXXXXXXXXXX**XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX**XXXXXXXXHFL1391 

ENTRY  FLUXE2(L,I,MIN)  HFL1392 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxHFL1393 

FI1=R0X*UX  HFL1394 

FI2=R0X*UXXX2+PX  HFL1395 

FI2=FI2-PX  HFL1396 

FI3=UXX(G12*PX+0.5XROX*UXXX2+ZXXROXXQDET)  HFL 1397 

FI4=ZXXR0XXUX  HFL1398 

ROU00=ROX*UX  HFL1399 

GO  T0(10,20,30,40,50,60),  NFLUX  HFL1400 

10  CONTINUE  HFL1401 
60    CONTINUE  HFL1402 

DFDXI1=DRDXIXXUX+R0XXDUDXIX  HFL 1403 

DFDXI2=DRDXIXXUXXX2+2.XR0X*UXXDUDXIX+DPDXIX  HFL 1404 

DFDXI2=DFDXI2-DPDXIX  HFL1405 

DFDXI3=DUDXIXX(G12XPX+0.5XROXXUXXX2)  HFL 14  06 

1      +UXX(G12XDPDXIX+0.5*DRDXIXXUXXX2+ROX*UXXDUDXIX)  HFL1407 

DFIDT1=DR0DTXXUX+R0XXDUDTX  HFL 14 08 

DFIDT2=DR0DTX*UXX*2+2.XR0XXUXXDUDTX+DPDTX  HFL 14 09 

DFIDT2=DFIDT2-DPDTX  HFL1410 

DFIDT3=DUDTXX(G12xPX+0.5xROXXUXXX2)  HFL 1411 

1      +UXX(G12XDPDTX+0.5*DRODTXXUXX*2+ROX*UXXDUDTX)  HFL1412 

FIDOTl=-ROU0  0xDFDXIl+DFIDTl  HFL 141 3 

FIDOT2=-ROU0  0*DFDXI2+DFIDT2  HFL 141 4 

FIDOT3=-ROU0  0XDFDXI3+DFIDT3  HFL 141 5 

UXDOT=-ROU0  0*DUDXIX+DUDTX  HFL 1416 

PXDOT=-ROUOOXDPDXIX+DPDTX  HFL 1417 

GXDOT=-ROU0  0XDGDXIX+DGDTX  HFL 1418 

ROXDOT=-ROU0  0*DRDXIX+DRODTX  HFL 1419 

FIH1=FI1+DT2XFID0T1  HFL1420 

FIH2=FI2+DT2XFID0T2  HFL1421 

GIH=PX+DT2XPXD0T  HFL1422 

FIH3=FI3+DT2XFID0T3  HFL1423 

FIH4=FI4  HFL1424 

UXN=UX+DTXUXDOT  HFL1425 

PXN=PX+DT*PXD0T  HFL1426 

GXN=GX+DTXGXD0T  HFL1427 

R0XN=R0X+DTXR0XD0T  HFL1428 

GO  TO  90  HFL1429 

20    CONTINUE  HFL1430 

EVO=GL*DSQRT(BETAO)  HFL1431 

201  CONTINUE  HFL1432 

DFIDA1=DR0DAX*UX+R0XXDUDAX  HFL 1433 

DFIDA2=DR0DAXXUXXX2+2.XR0X*UXXDUDAX+DPDAX  HFL 1434 

DFIDA2=DFIDA2-DPDAX  HFL1435 

DFIDA3=DUDAXX(G12XPX+0.5XROXXUXXX2)  HFL 1436 

1       +UXX(G12*DPDAX+0.5*DRODAXXUXXX2+ROXXUX*DUDAX)  HFL1437 

FID0T1=-EV0XDFIDA1  HFL1438 

FIDOT2=-EV0*DFIDA2  HFL1439 

FID0T3=-EV0XDFIDA3  HFL1440 
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E:  HFL       SCRIPT    Al 

FIH1=FI1+DT2*FID0T1  HFL1441 

FIH2=FI2+DT2*FID0T2  HFL1442 

FIH3=FI3+DT2XFID0T3  HFL1443 

FIH4=FI4  HFL1444 

GA=DGDAX  HFL1445 

IF(NFLUX.EQ.5)GA=-GA  HFL1446 

DROUA=UX*DRODAX+ROX*DUDAX  HFL1447 

BETAPR=0.5*DSQRT(BETA0)x(GA-DROUA)  HFL1448 

FIH2=FIH2-DPDBXXBETAPR*DT2  HFL1449 

UXDOT=-EV0*DUDAX+BETAPR*DUDBX  HFL 1450 

PXDOT=-EV0*DPDAX+BETAPR*DPDBX  HFL 1451 

GIH=PX+DT2*PXD0T  HFL1452 

GXDOT=-EV0*DGDAX+BETAPRXDGDBX  HFL 1453 

ROXDOT=-EV0*DRODAX+BETAPRXDRODBX  HFL 1454 

UXN=UX+DT*UXDOT  HFL1455 

PXN=PX+DT*PXDOT  HFL1456 

GXN=GX+DTXGXDOT  HFL1457 

ROXN=ROX+DT*ROXDOT  HFL1458 

GO  TO  90  HFL1459 

50    CONTINUE  HFL1460 

EVO=-GR*DSQRT(BETAO)  HFL1461 

GO  TO  201  HFL1462 

30    CONTINUE  HFL1463 

40    CONTINUE  HFL1464 

GO  TO  60  HFL1465 

90    CONTINUE  HFL1466 

RETURN  HFL1467 

END  HFL1468 
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